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Abstract 

This  report  documents  a  FORTRAN  version  of  the  Spencer-Mpiier-Weore  chemicoi 
thermodynamic  modei  for  aqueous  eiectrolyte  soiutions  at  subzero  temperatures 
(FREZCHEM).  FREZCFIEM  is  structured  to  predict  the  chemicoi  composition  and 
unfrozen  water  of  aqueous  soiutions  between  -60  °C  and  +25  °C  at  atmospheric 
pressure  (0.101325  MPa).  FREZCHEM  inciudes  two  reaction  pathways:  1) 
freezing  at  variabie  temperature  and  fixed  total  water  and  2)  evaporation  at 
variable  water  and  fixed  temperature.  Activity  coefficients  and  the  activity  of 
water  are  calculated  using  the  Pitzer  equations,  which  are  valid  to  high  solution 
ionic  strengths  (=20  mol  kg-i).  Fifteen  chloride  and  sulfate  salts  of  sodium, 
potassium,  calcium,  and  magnesium  are  included  in  the  model.  Predicted  and 
experimental  measurements  of  solute  molalities  and  the  unfrozen  water  fraction 
during  seawater  freezing  are  in  good  agreement.  At -50  °C,  0.3%  of  seawater 
remains  unfrozen  with  99.7%  of  Na  and  95.5  %  of  Cl  having  precipitated  into 
one  of  four  salts.  FREZCHEM  should  find  many  applications  in  physicochemical 
studies  of  aqueous  solutions  and  freezing. 
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FREZCHEM:  A  Chemical-Thermodynamic  Model  for 
Aqueous  Solutions  at  Subzero  Temperatures 

GILES  M.  MARION  AND  STEVEN  A.  GRANT 


INTRODUCTION 

The  presence  of  solutes  lowers  the  freezing- 
point  of  aqueous  solutions;  and  when  solutions 
freeze,  ice  formation  largely  excludes  solutes,  con¬ 
centrating  them  in  the  tmfrozen  brine  which,  in 
turn,  further  lowers  the  brine  freezing  point.  For 
example,  when  the  temperature  of  a  1.0  molal 
NaCl  solution  is  lowered  below  0  °C,  ice  starts  to 
form  at  approximately -3  °C  (Fig.  1).  As  ice  forms, 
the  NaCl  molality  increases  and  the  temperature 
drops  along  the  ice-solution  line,  assuming  that 
equilibrium  is  maintained,  imtil  the  eutectic  is 
reached  at  -21.3  °C  and  5.17  molal  (Fig.  1).  At  this 
point,  the  residual  solution  freezes,  forming  a 
mixture  of  ice  and  hydrohalite  [NaCl«2H20(cr)]. 
Clearly,  complex  interactions  exist  between  solute 
concentrations  and  the  freezing  process. 

A  number  of  chemical  thermodynamic  models 
for  aqueous  electrolyte  solutions  have  been  devel¬ 
oped  in  recent  years  (Sposito  and  Mattigod  1979, 
Plummer  et  al.  1988,  Spencer  et  al.  1990).  How¬ 
ever,  only  the  Spencer-MoUer-Weare  model  deals 
explicitly  with  aqueous  solutions  at  subzero  tem¬ 
peratures.  This  model  accounts  for  the  precipita¬ 
tion  of  the  chloride  and  sulfate  salts  of  sodium, 
potassium,  calcium,  and  magnesium  over  the  tem¬ 
perature  range  from  -60  °C  to  +25  °C  at  atmo¬ 
spheric  pressure  (0.101325  MPa).  While  the  pa¬ 
rameters  used  in  the  model  have  been  published 
(Spencer  et  al.  1990),  a  working  version  of  the 
model  has  not  been  published.  Such  a  working 
model  would  be  extremely  useful  for  geochemists 
and  geophysicists  interested  in  the  complex  inter¬ 
actions  between  solutes  and  the  freezing  process. 
This  preliminary  report  will  document  a  FOR¬ 
TRAN  version  of  the  Spencer-MoUer-Weare  model, 
which  we  have  chosen  to  call  FREZCHEM. 


MODEL  STRUCTURE 

FREZCHEM  is  a  chemical  thermodynamic 
"equilibrium"  model.  This  model  will  calculate 
the  equilibrium  composition  of  aqueous  solutions 
at  specified  temperatures,  but  will  not  provide  any 
information  on  the  time  required  to  reach,  or  the 
path  of  change  to,  the  equilibrium  state. 
FREZCHEM  calculates  equilibrium  between  wa¬ 
ter  and  ice  and  between  dissolved  and  solid  phase 
salts. 

For  water  and  ice  at  equilibrium, 

hi  -  l'^®H20(cr,I))  —  hw  =  hw°  + 

RT  ln(flH20(l))  (1) 

where  R  is  the  gas  constant,  T  is  absolute  tempera¬ 
ture,  Pi  and  Pw  are  the  chemical  potentials  of  ice 
and  water,  pi°  and  Pw°  are  the  standard  chemical 
potentials,  and  flH20(crl)  and  aH20(l)  are  the  activi¬ 
ties  of  water  in  the  form  of  ice  and  liquid  water. 
Equation  1  can  be  rearranged  to  yield 

=  In  {«H20(l)/«H20(cr,I))  = 

RT  '  RT 

(2) 

where  AfusG°  is  the  standard  Gibbs  energy  of 
fusion  (pi°  -  Pw°)-  At  equilibrium  with  respect  to 
pure  ice,  flH20(cr,  I)  =  l-O-  Therefore 

(flH20(l))  =  exp[AfosG7RT]  =  K  (3) 

the  activity  of  water  is  equal  to  the  equilibrium 
constant  (K)  at  a  given  temperature.  For  equilib¬ 
rium  with  respect  to  a  salt  (e.g.,  hydrohalite) 
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Figure  1.  A  stability  diagram  for  the  NaCl-H20  system  at  subzero  temperatures. 
Experimental  measurements  are  from  Bukshtein  et  al.  (1953)  and  Hall  et  al.  (1988). 


NaCl*2H20(cr)  ^ 

Na+(aq)  +  Cl-(aq)  +  2  H20(l)  (4) 

and 

(«Na+(aq))  {«a-(aq))  («H20(1))^  =  ^  (5) 

assuming  that  the  activity  of  the  pure  salt 
[NaCl»2H20(cr)  =  1.0.  The  relation  between  activ¬ 
ity  (a)  and  molality  (m)  of  a  substance  "B"  is  given 
by 

«B  =  yb  m  (6) 

where  y  is  the  activity  coefficient.  If  molalities 
are  given,  to  calculate  activities  requires  an 
extra-thermodynamic  model  for  activity  coeffi¬ 
cients. 

The  Pitzer  equations  are  used  to  calculate  the 
activity  of  water  and  the  activity  coefficients  for 
soluble  ions.  These  equations  allow  calculation  of 
activity  coefficients  to  high  ionic  strengths  in  com¬ 
plex  solutions  (Plummer  et  al.  1988,  Spencer  et  al. 
1990,  Pitzer  1991). 

The  osmotic  coefficient  ((})),  and  single-ion  activ¬ 
ity  coefficients  for  cations  (ym)  and  anions  (yx)  are 
given  by 


((j)  -  1)  =  2/(Znii)j-Al>Zl  5/(l  +  b  /O  S)  + 

IZm(Wa(6'l’ca  +  2:Cca) 

+  +  Sma  To-'a)  + 

ZZmam,(o'^aa'+2mc'Paa'c))  (7) 

In  (ym)  =  f  +  2:ma{2BMa  +  ZCMa)  + 

+  2^Wa^Mca) 

+ZZmama-'PaaM  +  |ziv^  SZOTcmaQa  (8) 

In  (yx)  =  Zx^  F  +  +  ZC^  + 

Zma(2®xa  +  ^mf¥xac) 

+ZZm  (W  c''**  ccX  +  |zx|  (W  aQa  (9) 

where  A*!*  is  the  Debye-Huckel  constant,  b  is  a 
constant  (=1 .2),  B't>,  C,  $'!>,  and  'P  are  salt  interaction 
coefficients,  z,  is  the  charge  number  of  the  ith  ion, 
I  is  the  ionic  strength  defined  by 

I  -  0.5 'Em  ,  (10) 

Z  is  defined  by 
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Table  1.  Minerals  currently  in  the  FREZCHEM 
model. 


1.  Na2SO4«10H2O(cr)  (mirabilite) 

2.  Na2S04(cr)  (thenardite) 

3.  K2S04(cr)  (arcanite) 

4.  Mg;S04*6H20(cr)  (hexahydrite) 

5.  MgS04*7H20(cr)  (epsomite) 

6.  MgS04*K2S04  6H20(cr)  (picromerite) 

7.  NaCl(cr)  (halite) 

8.  NaCl*2H20(cr)  (hydrohalite) 

9.  KCt(cr)  (sylvite) 

10.  CaCl2*6H20(cr)  (antarcticite) 

11.  MgCl2»6H20(cr)  (bischofite) 

12.  MgCl2*8H20(cr) 

13.  MgCl2-12H20(cr) 

14.  KMgCl3*6H20(cr)  (camallite) 

15.  CaCl2*2MgCl2*12H20(cr)  (tachyhydrite) 


Z  =  -Lm  i\zi\  (11) 

and  F  is  a  complex  function  of  A^,  I,  and  the  salt 
interaction  coefficients  (Plummer  et  al.  1988,  Pitzer 
1991).  Given  ^  from  eq  7,  the  activity  of  ivater  is 


defined  by 

flH20(l)  =  exp(-  (j)  Zmi/55.50837) .  (12) 

The  equilibrium  constants  and  the  Pitzer  equa¬ 
tion  parameters  for  the  temperature  range  -60  °C 
to  +25  °C  were  published  by  Spencer  et  al.  (1990)  as 
functions  of  temperature  using  equations  of  the 
form 

P(T)  =  fli  +  02  r  +  fl6  +  «9  + 

a^/T  +  04  ln(T)  (13) 


where  P(T)  is  the  Pitzer  equation  parameter  di¬ 
rectly  or  the  function,  [-ArG°/RT]  (where  Ai.G°  is 
the  standard  reaction  Gibbs  energy),  for  equilib¬ 
rium  constants  [fC  =  exp(-ArG°  / J?T)]  and  T  is  ther¬ 
modynamic  temperature  (K).  Refer  to  Spencer  et 
al.  (1990)  for  the  specific  parameters  used  in  the 
model  (their  Tables  1,  2,  and  3).  The  only  addi¬ 
tional  parameter,  beyond  those  used  in  the  Spen- 
cer-MoUer-Weare  model,  is  T(Ca,Mg,S04) =0.024, 
which  is  treated  as  a  constant  independent  of 
temperature  (Pitzer  1991). 


Table  2.  A  listing  of  chemical  species  and  their  numerical  designation  in  the  FREZCHEM  model. 


A.  Solution  and  atmospheric  species 


# 

Species 

# 

1 

Na+(aq) 

11 

2 

K+(aq) 

12 

3 

Ca2+(aq) 

13 

4 

Mg2+(aq) 

14 

5 

H+(aq) 

15 

6 

16 

7 

17 

8 

18 

9 

19 

10 

20 

Species 

# 

Species 

Cl-(aq) 

21 

C02(aq) 

S042-(aq) 

22 

CaS04°(aq) 

OH-(aq) 

23 

MgS04‘'(aq) 

HC03-(aq) 

24 

COsMaq) 

25 

26 

27 

28 

29 

C02(g) 

30 

H20(1) 

B.  Solid  phase  species 


# 

Species 

# 

Species 

it 

31 

H20(cr,l) 

41 

Na2SO4*10H2O(cr) 

51 

32 

NaCl-2H20(cr) 

42 

Na2S04(cr) 

52 

33 

NaCl(cr) 

43 

MgS04*6H20(cr) 

53 

34 

KCl(cr) 

44 

MgS04  •7H20(cr) 

54 

35 

CaCl2«6H20(cr) 

45 

K2S04(cr) 

55 

36 

MgCl2*6H20(cr) 

46 

MgS04»K2S04*6H20(cr) 

56 

37 

MgCl2'8H20(cr) 

47 

57 

38 

MgCl2*12H20(cr) 

48 

58 

39 

KMgCla  •6H20(cr) 

49 

59 

40 

CaCl2*2MgCl2‘12H20(cr) 

50 

60 
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Table  3.  FREZCHEM  model  output  for  freezing  of  seawater  at  -30  °C. 


Temp  (K) 

243.15 

Ion.  Sir. 

7.8770 

AH20 

0.74696 

Phi 

1.6390 

H2O  (g) 

33.120 

Ice  (g) 

948.20 

Solution 

Initial 

Final 

Species 

Cone. 

Cone. 

Aet.  Coef. 

Aetivity 

Moles 

Na 

0.48695 

1.6940 

0.54243 

0.91885 

0.56104E-01 

K 

0.10630E-4)1 

0.32095 

0.19599 

0.62905E-01 

0.10630E-01 

Ca 

0.95300E-02 

0.28673 

0.51732 

0.14833 

0.94967E-02 

Mg 

0.55160E-01 

1.6654 

1.1282 

1.8790 

0.55160E-01 

Cl 

0.56818 

5.9081 

1.8095 

10.691 

0.19568 

SO4 

0.29390E-01 

0.55838E-02 

0.32328 

0.18051E-02 

0.18494E-03 

CaS04 

0.00000 

0.10058E-02 

1.0000 

0.10058E-02 

0.33313E-04 

MgS04 

0.00000 

0.42492E-05 

1.0000 

0.42492E-O5 

0.14074E-06 

HzOd) 

0.74696 

1.8384 

Solid 

Equil. 

Species 

Moles 

Constant 

Ice 

52.633 

0.74693 

NaCl2-2H20 

0.37250 

5.4770 

NaCl 

0.00000 

19.127 

KCl 

0.00000 

1.0399 

CaCl2*6H20 

0.00000 

695.12 

MgCl2-6H20 

0.00000 

40585.0 

MgCl2*8H20 

0.00000 

1403.3 

MgCl2*12H20 

0.00000 

27.775 

KMgCl3*6H20 

0.00000 

1054.2 

CaCl2*2MgCl2*12H20 

0.00000 

0.17775E+21 

Na2SO4-10H2O 

0.29172E-01 

0.81213E-04 

Na2S04 

0.00000 

0.73236 

MgS04  *61120 

0.00000 

3.1211 

MgS04  *71120 

0.00000 

0.93967E-O1 

K2SO4 

0.00000 

0.61722E-02 

MgS04*K2S04*6H20 

0.00000 

0.21935E-O4 

CaS04  (ion-pair) 

0.26618 

MgS04  (ion-pair) 

798.18 

Iterations  =  39 


The  specific  salts  included  in  the  present  model 
are  listed  in  Table  1.  In  addition  to  these  salts,  the 
model  also  explicitly  considers  ion-pairs  for  Ca- 
SO4  and  MgS04. 

MATHEMATICAL  ALGORITHM 

Two  basic  mathematical  approaches  are  used  in 
chemical  thermodynamic  models  to  calculate  equi¬ 
librium  concentrations  and  activities:  Gibbs  en¬ 
ergy  minimization  or  solution  of  sets  of  nonlinear 
equations.  The  Spencer  et  al.  (1990)  model  explic¬ 
itly  uses  Gibbs  energy  minimization.  The 


PHRQPITZ  model  (Plummer  et  al.  1988)  uses  the 
nonlinear  equation  solution  method.  Both  ap¬ 
proaches  ultimately  provide  the  same  answers. 

The  FREZCHEM  model  solves  sets  of  nonlinear 
equations  that  include  ice-water  and  mineral 
equilibria.  These  equations  may,  in  principle,  be 
solved  either  simultaneously  or  sequentially. 
FREZCHEM  uses  the  sequential  approach,  which 
treats  each  reaction  as  an  independent  reaction 
(Fig.  2).  The  solution  of  the  set  of  nonlinear  equa¬ 
tions  occurs  by  iterating  through  the  sequence 
many  times.  Advantages  of  the  sequential  ap¬ 
proach  vis-a-vis  the  simultaneous  or  Gibbs  energy 
minimization  approaches  include  1)  programming 
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Figure  2.  The  FREZCHEM  model  flowchart. 


simplicity  since  only  one  reaction  is  considered  in 
separate  subroutines,  2)  ease  of  adding  or  sub¬ 
tracting  new  reactions  since  they  can  be  inserted  or 
deleted  at  any  point  in  the  sequence,  and  3)  the 
ease  of  program  operation  for  both  simple  and 
complex  aqueous  solutions.  With  respect  to  the 
last  point,  only  subroutines  containing  the  ele¬ 
ments  in  the  "input"  listing  are  used.  For  example, 
if  only  Na  and  Cl  are  present,  then  the  program 
wiU  only  "call"  the  subroutine  dealing  with  NaCl 
equilibrium  and  will  bypass  the  remainder.  For 
complex  aqueous  solutions  such  as  seawater,  the 
program  will  "call"  every  subroutine.  The  major 
disadvantage  of  the  sequential  approach  is  that  it 
can  be  slow  to  converge.  Another  problem  with 
the  current  model  is  convergence.  At  times,  the 
program  diverges  or  oscillates  and  never  arrives  at 
the  equilibrium  concentrations. 

The  program  assumes  that  the  initial  aqueous 
solution  is  imfrozen  and  calculates  the  amount  of 
ice,  if  any,  that  will  form  at  the  specified  salt  con¬ 


centrations  and  temperature.  The  program  first 
tests  to  see  if  ice  and  water  are  at  equilibrium  at  the 
specified  temperature.  The  criteria  are 

IF  flH20(l)  >  ^  THEN  freeze  (14) 

IF  ^JH20(1)  =  ^  THEN  equilibrium  (15) 

IF  «H20(1)  <  ^  THEN  water  is  stable  phase 

(16) 

where  aH20(l)  is  the  activity  of  water  calculated 
from  the  Pitzer  equations  (eq  7  and  12)  and  K  is  the 
equilibrium  constant  for  ice-water  (eq3).  If  flH20(l) 
>  K,  this  implies  that  water  is  thermodynamically 
rmstable  and  ice  (stable  phase)  should  form.  Ice 
formation  is  calculated  indirectly.  First,  solute  con¬ 
centrations  are  increased  proportionally,  which 
decreases  flH20(l)  (eq  12)  imtil  aH20(l)  =  arid 
equilibrium  is  reached.  Ice  formation  is  calculated 
by 
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Ice  =  TotalH20  -  Watetf  (17) 

Watetf  =  Water;  (Ewi/XjWf)  (18) 

where  Waterf,  Water;,  and  TotalHaO  are  the  final 
water,  initial  water,  and  total  H2O  (ice  and  water), 
m;  are  the  initial  molahties  of  solutes,  and  m{  are 
the  final  molahties.  These  calculations  assume  that 
the  ice  phase  is  a  pure  phase.  If  aH20(l)  <  K  this 
imphes  that  water  is  the  thermod3mamicaUy  stable 
phase  and  ice  should  not  form  at  the  specified  salt 
concentrations  and  temperature. 

In  addition  to  freezing,  FREZCHEM  will  also 
evaporate  solutions.  In  contrast  to  freezing,  the 
evaporation  of  water  is  not  controUed  by  an  equi¬ 
librium  relation.  The  evaporation  routine  simply 
removes  water  incrementally  from  the  solution,  at 
a  fixed  temperature,  according  to  the  user's  speci¬ 
fications.  Tlie  program  will  not  allow  evaporation 
if  ice  is  thermod3mamically  stable  at  the  specified 
temperature  and  solution  concentrations. 

After  ice  and  water  are  equhibrated,  the  pro¬ 
gram  cycles  through  subroutines  to  handle  min¬ 
eral  equilibria  (Fig.  2).  At  present,  there  are  11 
subroutines  to  handle  the  15  mineral  phases  (Table 
1)  and  the  two  ion-pair  equihbria.  Minerals  that 
are  identical  in  chemical  composition  except  for 
waters  of  hydration  are  handled  in  the  same  sub- 
rbutme  [e.g.,  NaCl(cr)  and  NaCl«2H20(cr)]. 

After  cycling  through  the  mineral  equilibrium 
equations,  the  program  calculates  the  amount  of 
water  associated  with  hydrated  salts  (Fig.  2) .  Then 
the  program  tests  for  overall  mineral  equilibri¬ 
um.  In  FREZCHEM,  the  initial  and  final  solution 
phase  concentrations  must  agree  to  within  ±0.1  % 
for  all  solution  phase  cations  and  anions.  In  gen¬ 
eral  if  more  than  one  mineral  phase  is  precipitat¬ 
ing,  the  first  cycle  through  the  mineral  subrou¬ 
tines  will  not  suffice  to  establish  overall  mineral 
equilibrium  and  the  equations  must  be  solved 
repeatedly. 

Eventually  mineral  equilibrium  is  established. 
Then  the  program  tests  to  see  if  ice  and  water  are 
stiU  in  equilibrium  or  aH20(l)  <  K.  For  equilibrium, 
aH20(l)  must  agree  with  K  to  within  ±0.00005.  If  ice 
and  water  are  not  in  equilibrium  and  aH20(l)  > 
then  the  program  cycles  through  the  ice-water 
equilibrium  and  subsequently  mineral  equilib¬ 
rium  subroutines.  Eventually  the  program  con¬ 
verges  to  within  the  numerical  equilibrium  crite¬ 
ria  established  in  the  program,  and  output  is  sent 
to  the  print  files. 


FREZCHEM  PROGRAM 

A  listing  of  the  FREZCHEM  FORTRAN  pro¬ 
gram  is  in  Appendix  A.  A  major  effort  was  made 
to  adhere  to  the  FORTRAN  77  ANSI  Standard  in 
writing  this  program.  For  example,  FREZCHEM 
uses  an  80-column  statement  line  where  columns 
1-5  are  statement  labels,  column  6  is  a  continua¬ 
tion  column,  columns  7-72  are  used  for  state¬ 
ments,  and  columns  73-80  are  used  for  commen¬ 
tary.  Variable  names  are  <  6  characters. 
FREZCHEM  is  a  FORTRAN  translation  of  an  ear¬ 
lier  program  written  in  TrueBASIC.  FORTRAN 
assumes,  unless  otherwise  declared,  that  variables 
beginning  with  the  letters  "IJKLMN"  are  integer 
variables  and  all  other  initial  letters  refer  to  real 
variables.  In  order  to  retain  the  same  variable 
names  between  programs,  this  restriction  required 
formal  declaration  of  some  real  variables  in  the 
FORTRAN  version.  FREZCHEM  is  currently  nm- 
ning  on  a  Macintosh  Quadra  800  computer  using 
the  MacFortran  II  compiler  (version  3.2,  Absoft 
1991). 

The  FREZCHEM  program  consists  of  a  main 
program  called  FREZCHEM  and  15  subroutines. 
FREZCHEM  reads  in  data  from  screen  queries, 
initializes  some  parameters,  freezes  or  evaporates 
the  aqueous  solution,  calculates  the  water  of  hy¬ 
dration,  tests  for  chemical  thermodynamic  equi¬ 
libria,  and  most  importantly  "calls"  the  various 
subroutines  (Fig.  2). 

Subroutine  PARAMETER  calculates  the  param¬ 
eters  for  the  Pitzer  equations  and  the  equilibrium 
constants  as  functions  of  temperature  (eq  13).  For 
the  most  part,  these  parameters  are  taken  from 
Spencer  et  al.  (1990).  Subroutine  PITZER  calcu¬ 
lates  activity  coefficients  and  the  activity  of  water 
using  the  Pitzer  and  other  equations  (eq  7-12). 
Subroutine  INTERACTION  calculates  the  higher- 
order  electrostatic  interaction  terms  for  the  Pitzer 
equations.  Subroutine  PPRINT  prints  the  results. 
The  remaining  11  subroutines  deal  with  mineral 
equilibria .  For  example.  Subroutine  MGCL2  deter¬ 
mines  which  MgCl2  salt  is  thermod3mamically 
stable  and  equilibrates,  if  necessary,  the  salt  with 
the  solution  phase.  Other  mineral  subroutines 
perform  similar  functions  for  specific  minerals. 

FREZCHEM  was  structured  to  be  user-expand¬ 
able.  To  alter  the  model  requires  knowing  how 
chemical  species  are  expHcitly  designated  in  the 
model.  Table  2  summarizes  the  chemical  species 
currently  in  or  anticipated  to  be  included  in  the 
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model  in  the  near  future.  Species  1-10  are  reserved 
for  soluble  cations;  species  11-20  are  reserved  for 
soluble  anions;  species  21-30  are  reserved  for  neu¬ 
tral  species,  both  solution  phase  and  gas  phase. 
The  solid  phase  species  are  designated  with  num¬ 
bers  31-60.  Blank  spaces  in  Table  2  are  available  for 
future  additions. 


PROGRAM  INPUT  AND  OUTPUT 

Input  to  FREZCHEM  is  through  screen  queries 
that  request  the  molal  concentrations  for  Na,  K, 
Ca,  Mg,  Cl,  and  SO4.  FREZCHEM  is  structured  to 
cycle  through  a  sequence  of  temperature  changes 
(freezing)  or  water  removals  (evaporation).  For  a 
temperature  change,  the  screen  queries  request 
the  initial  temperature,  the  final  temperature,  and 
the  incremental  temperature  change.  For  example, 
to  calculate  the  freezing  of  seawater  between 273.15 
K  (0  °C)  and  223.15  K  (-50  °C)  by  5  K  intervals, 
enter:  273.15,  223.15,  and  5.0.  Because  the  model 
assumes  an  initially  unfrozen  solution,  enter  the 
temperature  sequence  from  high  (initial)  to  low 
(final).  For  evaporative  removal,  the  screen  que¬ 
ries  request  the  temperature  (assumed  fixed  for 
evaporation),  the  final  total  water  content  (grams), 
and  the  incremental  water  removal  (grams).  To 
evaporate  water  from  1000  g  (the  fixed  initial 
water  content  for  both  freezing  and  evaporation) 
to  100  g  by  50-g  increments  at  0.0  °C,  enter  273.15, 
100,  and  50  for  the  evaporative  screen  queries. 

Output  from  the  program  is  to  a  file  called 
"FrData"  which  can  be  sent  to  a  printer.  A  freezing 
simulation  for  seawater  at  243.15  K  (-30  °C)  is 
given  in  Table  3,  and  an  evaporative  (1000  g  to  100 
g)  simulation  for  seawater  at  0  °C  is  given  in  Table 
4.  After  the  title,  the  first  line  of  output  gives  the 
temperature,  ionic  strength,  activity  of  water 
(AH2O),  osmotic  coefficient  (phi),  and  the  final 
calculated  amounts  of  water  and  ice.  The  basic 
tmit  of  mass  in  the  model  is  1.0  kg  water.  For  the 
seawater  simulation  at  -30  °C,  the  final  distribu¬ 
tion  of  water  is  948.20  g  of  ice,  33.12  g  water,  13.42 
g  of  water  in  hydrohalite  (0.3725  moles),  and  5.26 
g  of  water  in  mirabilite  (0.0292  moles),  which  adds 
up  to  exactly  1000.00  g  (T able  3) .  The  final  distribu¬ 
tion  of  water  in  the  evaporative  scenario  at  0  °C  is 
96.93gwaterand3.07gwaterin  mirabilite  (0.017043 
moles),  which  adds  up  to  exactly  100.00  g  (Table  4). 

The  solution  phase  species  are  grouped  together 
in  the  output.  This  output  includes  initial  concen¬ 


tration  (the  program  input),  final  concentration, 
activity  coefficient,  activity,  and  the  moles  remain¬ 
ing  in  the  solution  phase.  Because  the  basic  imit  of 
mass  is  1.0  kg  water,  the  initial  concentration 
[moles/kg  (water)]  is  also  equal  to  the  total  moles. 
A  difference  between  "Initial  Cone."  and  "Moles" 
in  the  solution  phase  output  is  attributable  to 
either  salt  precipitation  or  ion-pair  formation.  For 
example,  Ihe  initial  concentration  (total  moles)  of 
Na  is  0.48695.  The  distribution  of  Na  at  -30  °C  is 
0.05610  moles  in  the  solution,  0.37250  moles  as 
hydrohaHte,  and 0.05834 moles  as  mirabilite,  which 
adds  up  to  0.48694  moles  of  Na  (Table  3).  For  the 
seawater  evaporative  simulation,  the  final  distri¬ 
bution  of  Na  is  0.45286  moles  in  the  solution  and 
0.03409  moles  as  mirabilite,  which  adds  up  to 
0.48695  moles  of  Na  (Table  4) 

The  solid  phase  species  are  also  grouped  to¬ 
gether  in  the  output.  This  output  includes  the 
moles  of  solid  phase  species  and  the  equilibrium 
constant  for  the  reaction.  For  the  seawater  freezing 
simulation,  ice,  hydrohalite,  and  mirabilite  have 
precipitated  at  -30  °C  (Table  3).  For  the  seawater 
evaporative  simulation,  only  mirabilite  appears  to 
have  precipitated  at  0  °C  (Table  4);  however,  see 
Program  Limitations  (below)for  a  discussion  of  prob¬ 
able  gypsum  precipitation.  Whenever  ice  is  present, 
the  activity  of  water  (0.74696)  must  agree  with  the 
equilibrimn  constant  for  ice-water  (0.74693)  (Table 
3)  to  within  ±0.00005  units,  which  is  the  conver¬ 
gence  criterion  used  in  the  model.  For  salts  that  are 
precipitating,  the  appropriate  solution  phase  ac¬ 
tivity  product  should  equal  the  equilibrium  con¬ 
stant.  For  NaCl  •2H20(cr),  the  activity  product  is 
equal  to  (0.91885)  x  (10.691)  x  (0.74696)2  =  5.48 
which  agrees  with  the  equilibrium  constant  of  5 .48 
(Table  3). 

The  final  output  are  the  equilibrium  constants 
for  the  CaS04  and  MgS04  ion-pairs  and  the  num¬ 
ber  of  iterations  through  the  equilibrium  routines 
before  convergence. 

PROGRAM  LIMITATIONS 

FREZCHEM  has  convergence  problems  espe¬ 
cially  at  high  ionic  strengths  (>15  molal),  where 
activity  coefficients  are  rapidly  changing,  and  at 
junctions,  where  new  phases  begin  to  precipitate. 
The  program  contains  a  routine  to  catch  math¬ 
ematical  problems  or  special  circumstances.  If  "it¬ 
erations"  reaches  300,  a  message  to  this  effect  is 
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Table  4.  FREZCHEM  model  output  for  evaporation  of  seawater  at  0  °C. 


Temp  (K) 

273.15 

Ion.  Str. 

6.8871 

AH20 

0.75383 

Phi 

1.3728 

H2O  (g) 

96.930 

Ice  (g) 

0.00000 

Solution 

Initial 

Final 

Species 

Cone. 

Cone. 

Act.  Coef. 

Activity 

Moles 

Na 

0.48695 

4.6721 

0.81929 

3.8278 

0.45286 

K 

0.10630E-01 

0.10967 

0.39751 

0.43593E-01 

0.10630E-01 

Ca 

0.95300E-02 

0.92284E-01 

1.0610 

0.97912E-01 

0.89451E-02 

Mg 

0.55160E-01 

0.56905 

2.1926 

1.2477 

0.55158E-01 

Cl 

0.56818 

5.8618 

1.2301 

7.2105 

0.56818 

SO4 

0.29390E-01 

0.12133 

0.30951E-01 

0.37552E-02 

0.11760E-01 

CaS04 

0.00000 

0.60346E-4)2 

1.0000 

0.60346E-4)2 

0.58493E-03 

MgS04 

0.00000 

0.19557E-04 

1.0000 

0.19557E-04 

0.18957E-05 

HzCKl) 

0.75383 

5.3804 

Solid 

Equil. 

Species 

Moles 

Constant 

Ice 

0.00000 

1.0002 

NaCl-2H20 

0.00000 

17.957 

NaCl 

0.00000 

31.339 

KCl 

0.00000 

3.7731 

CaCl2*6H20 

0.00000 

1862.7 

MgCl2*6H20 

0.00000 

55958. 

MgCl2*8H20 

0.00000 

7080.2 

MgCl2*12H20 

0.00000 

708.26 

KMgCl2-6H20 

0.00000 

7819.3 

CaCl2-2MgCl2*12H20 

0.00000 

0.27626E+19 

Na2SO4-10H2O 

0.17043E-01 

0.32600E-02 

Na2S04 

0.00000 

0.69007 

MgS04*6H20 

0.00000 

0.18782 

MgS04«7H20 

0.00000 

0.54675E-01 

K2SO4 

0.00000 

0.69458E-D2 

MgS04*K2S04*6H20 

0.00000 

0.50713E-4)4 

CaS04  (ion-pair) 

0.60935E-01 

MgS04  (ion-pair) 

239.57 

Iterations  =  3 

printed  and  the  program  exits.  The  most  likely 
explanations  for  nonconvergence  within  300  itera¬ 
tions  are  divergence,  oscillating  around  the  cor¬ 
rect  solution,  or  the  system  is  beyond  the  eutectic, 
at  which  point  only  solid  phases  are  thermody¬ 
namically  stable. 

To  rectify  mathematical  problems,  try  changing 
the  temperature  or  evaporative  increments,  either 
increasing  or  decreasing  their  step  size  (e.g.,  change 
temperature  decrement  from  10  to  5  degrees  or 
change  water  decrement  from  50  to  100  g). 

Generally  one  can  tell  the  approximate  eutectic 
of  more  complex  mixtures  from  the  behavior  of 
pure  salts.  For  example,  the  eutectic  temperature 
for  pure  CaCl2  is  -50.4  °C  (Spencer  et  al.  1990). 
Failure  of  the  model  to  converge  at  ^0  °C  for  a 


solution  containing  CaCl2  such  as  seawater  carmot 
be  due  to  exceeding  the  eutectic,  which  must  be 
lower  in  temperature  than  the  eutectic  of  the  pure 
solution.  Failure  of  the  model  to  converge  within 
300  iterations  at  -54  °C  for  seawater  is  caused  by 
exceeding  the  eutectic  at  -53.6  °C  (Spencer  et  al. 
1990). 

If  you  use  FREZCHEM  to  identify  "junctions" 
such  as  the  eutectic,  you  need  to  carefully  scruti¬ 
nize  the  output  for  accuracy.  For  example, 
FREZCHEM  prints  results  for  seawater  freezing 
down  to  -53.7  °C;  at  -53.8  °C,  the  program  fails  to 
converge  within  300  iterations  because  the  solu¬ 
tion  is  beyond  the  eutectic  temperature.  If  one  ex¬ 
amines  the  printed  output  at  -53.7  °C,  the  solution 
is  slightly  supersaturated  with  respect  to  CaCl2 
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(the  last  of  the  seawater  salts  to  precipitate)  based 
on  calculated  activity  products,  but  the  program 
indicates  no  precipitation  of  CaCl2  (moles  =  0.0). 
At  -53.6  °C,  the  solution  is  imdersaturated  with 
respect  to  CaCl2.  The  true  eutectic  lies  between 
-53.7  °C  and  -53.6  °C,  which  agrees  with  Spencer 
et  al.  (1990)  who  place  the  seawater  eutectic  at 
-53.64  °C.  Small  errors  can  occur  at  junctions; 
therefore,  these  solutions  need  to  be  carefully  scru¬ 
tinized. 

FREZCHEM  is  limited  to  applications  where 
chloride  and  sulfate  salts  of  sodium,  potassium, 
calcium,  and  magnesium  are  dominant.  Two  com¬ 
mon  salts  not  currently  in  the  model  are  gypsum 
[CaS04»2H20(cr)]  and  calcite  [CaC03(cr)].  In  the 
evaporative  scenario  (Table  4),  the  calculated  ac¬ 
tivity  product  for  g3^sum  in  seawater  is  2.09  x  10^ 
which  is  apparently  supersaturated  with  respect 
to  gypsum  at  0°C{K  =  2.33  x  10“5,  Plummer  et  al. 
1988).  The  reason  gypsum  and  calcite  are  missing 
from  the  current  model  is  because  these  salts  are 
relatively  insoluble  and  as  a  consequence,  there 
are  few  solubility  data  at  subzero  temperatures. 
Because  neither  gypsum  nor  calcite  can  be  param¬ 
eterized  over  the  temperature  range  of  interest 
(-60  to  +25  °C),  we  have  chosen  not  to  include  their 
solubility,  even  over  a  limited  temperature  range. 


in  this  preliminary  version  of  FREZCHEM. 

This  is  not  to  say  that  FREZCHEM  cannot  be 
used  in  circumstances  where  gypsum  or  calcite  are 
precipitating,  only  that  caution  is  necessary  in 
interpreting  the  results.  For  example,  calcite  pre¬ 
cipitates  during  both  seawater  freezing  (Nelson 
and  Thompson  1954)  and  seawater  evaporation, 
as  does  gypsum  during  seawater  evaporation 
(McCaffrey  et  al.  1987).  However,  the  amoxmts  of 
calcite  precipitation  are  relatively  minor  and  can 
apparently  be  ignored.  Note  the  good  agreement 
in  Ca  concentrations  during  seawater  freezing 
between  experimental  data  and  the  model  that 
ignored  both  calcite  and  gypsum  (Figs.  3  and  4). 
One  cannot  ignore  gypsum  during  seawater  evapo¬ 
ration  (McCaffrey  et  al.  1987,  Herut  et  al.  1990) 
where  explicit  recognition  must  be  made  of  gyp¬ 
sum  precipitation.  As  pointed  out  previously,  a 
10-fold  concentration  of  seawater  at  0  °C  is  likely  to 
lead  to  the  significant  precipitation  of  both 
mirabUite  (0.0170  moles)  and  g5q)sum  (supersatu¬ 
rated)  (Table  4).  Where  gypsum  precipitation  is 
possible,  check  the  calculated  gypsum  activity 
product  for  apparent  supersaturation.  Unfortu¬ 
nately,  reliable  gypsum  solubility  products  are 
available  only  for  temperatures  >  0  °C  (Plummer  et 
al.  1988). 


Molality 


Figure  3.  The  concentrations  of  major  seawater  constituents 
during  freezing.  Experimental  measurements  are  from  Nelson 
and  Thompson  (1954). 
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freezing.  Experimental  measurements  are  from  Nelson  and  Thomp¬ 
son  (1954). 


PROGRAM  VALIDATION 

The  Spencer-M0ller-Weare  model  and  conse¬ 
quently  the  FREZCHEM  model  are  parameter¬ 
ized  primarily  with  data  from  pure  binary  and 
ternary  salt  solutiorrs.  Spencer  et  al.  (1990)  present 
13  figures  demonstrating  how  weU  their  model  fits 
the  experimental  data  (e.g..  Fig.  1).  In  general,  the 
agreement  is  excellent,  which  is  not  surprising 
because  this  is  the  database  used  to  parameterize 
the  model,  and  as  a  consequence,  this  is  not  a  true 
validation  of  the  model. 

The  primary  validation  of  the  Spencer-Moller- 
Weare  model  is  a  comparison  of  model  predic¬ 
tions  for  seawater  freezing  with  experimental  data 
from  Nelson  and  Thompson  (1954).  The  compari¬ 
son  of  major  (Fig.  3)  and  minor  (Fig.  4)  seawater 
constituents  demonstrates  good  agreement  except 
for  Mg  and  K  concentrations  between  -30°  and  -40 
°C  and  SO4  concentrations  between  -5°  and  -10 
°C.  The  seawater  experimental  data  between  -30° 
and  -40  °C  show  a  great  excess  of  cations  over 
anions  indicating  a  major  problem  with  the  ana¬ 
lytical  measurements.  All  of  the  discrepancy  be¬ 
tween  measured  and  predicted  Mg  and  K  concen¬ 
trations  can  be  attributed  to  this  charge  imbalance. 
The  model  predicts  that  mirabilite  precipitates  at 
-5.9  °C,  while  experimental  measurements  sug¬ 
gest  that  precipitation  begins  at  -8.2  °C.  This  is  the 


largest  discrepancy  in  the  temperature  of  first 
appearance  of  salt  between  the  model  and  the 
experimental  data.  We  plan  to  address  these  prob¬ 
lems  with  experimental  work  in  the  near  future. 

Another  approach  to  validating  the  model,  not 
used  by  Spencer  et  al.  (1990),  is  to  examine  the 
distribution  of  water  during  seawater  freezing. 
For  a  system  starting  with  1.0  kg  of  water  at  0  °C, 
approximately  50%  of  the  water  has  turned  to  ice 
by  the  time  the  temperature  has  dropped  to  -4  °C 
(Fig.  5).  By-10°C,  only  20%  of  the  original  water  is 
xmfrozen.  At  about  -23  °C,  hydrated  salts  become 
a  significant  sink  for  water;  this  is  where  hydrohalite 
begins  to  precipitate.  At  -50  °C,  the  original  1000  g 
of  water  is  redistributed  as  964.73  g  of  ice,  32.12  g 
of  hydrated  salts,  and  3.15  g  of  rmfrozen  water. 
There  is  good  agreement  between  calculated  un¬ 
frozen  water  and  experimental  measurements 
(Richardson  1976)  at  least  down  to  -36  °C  (Fig.  6). 
Below  -36  °C,  the  experimental  and  model  esti¬ 
mates  diverge  but  come  back  together  at  -50  °C 
where  the  model  predicts  3.15  g  water  and  the 
experimental  measurement  indicates  3.53  g  water. 
The  two  curves  diverge  at  the  point  where 
MgCl2*12H20(cr)  begins  to  precipitate  at  -36  °C. 
Given  major  rmcertainties  in  both  model  param¬ 
eterization  and  experimental  measurements,  we 
cannot  determine  if  the  model  or  experimental  mea¬ 
surements  are  more  accurate  at  low  temperature. 
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Figure  5.  The  model  prediction  of  the  distribution  of  water  during  seawater  freezing. 


The  importance  of  salt  precipitation  in  control¬ 
ling  seawater  freezing  is  illustrated  by  examining 
the  distribution  of  total  Na  and  Cl  during  the 
freezing  process  (Fig.  7).  At  -50  °C,  only  0.3  %  and 
4.5  %  of  the  original  soluble  Na  and  Cl  at  0  °C 
remain  in  the  solution  phase;  99.7  %  and  95.5  %  of 
the  Na  and  Cl  have  precipitated  in  four  salts. 


The  ionic  strength  of  the  imfrozen  seawater 
changes  from  0.72  to  11 .85  molal  between  0  °C  and 
-50  °C.  This  large  increase  in  salt  concentration 
during  the  freezing  process  is  the  primary  reason 
why  it  is  necessary  to  use  the  Pitzer  equations  to 
evaluate  the  activity  of  water  and  activity  coeffi¬ 
cients. 


ments  are  from  Richardson  (1976). 
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iii  Solution  H  NagSO^ 


Temperature  (°C) 

Figure  7.  The  model  predictions  of  sodium  and  chloride  distribution  during 
seawater  freezing. 


FUTURE  DIRECTIONS 

The  data  presented  in  this  report  (Fig.  3-7) 
demonstrate  the  great  complexity  of  the  solute¬ 
freezing  interaction  and  why  this  interaction  must 
be  considered  to  define  accurately  the  freezing 
process,  especially  at  high  salinities.  The 


FREZCHEM  model  can  calculate  the  role  of  the  Cl 
and  SO4  salts  of  Na,  K,  Ca,  and  Mg  in  the  freezing 
process.  The  model  is  useful  as  is,  but  more  work 
needs  to  be  done  in  order  to  develop  a  new  "state- 
of-the-art"  model. 

For  the  chloride  salts,  KCl  solubility  is  the  major 
imcertainty.  The  parameterization  of  KCl  solubil- 
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ity  by  Spencer  et  al.  (1990)  is  based  on  only  one  data 
point  below -24  °C.  We  are  beginning  experiments 
to  measure  the  solubility  of  KCl  in  KCl-CaCl2 
solutions  to  -50  °C.  At  a  minimum,  we  want  to  add 
gypsum  [CaS04»2H20(cr)]  solubility  to  the  model 
and  re-examine  the  solubility  of  Na2S04,  which  is 
not  predicted  well  by  the  current  model  (Fig.  4). 
And  finally,  we  need  to  add  calcite  [CaC03(cr)] 
solubility  to  the  model.  Calcite  and  gypsum  are 
common  minerals  that  are  present  in  many  natural 
systems.  Their  inclusion  in  the  model  would  be 
critical  to  many  applications. 

There  are  mathematical  convergence  problems, 
at  times,  with  the  FREZCHEM  model,  probably 
due  to  the  sequential  approach  to  solving  a  set  of 
nonlinear  equations  used  in  the  model.  An  alter¬ 
native  mathematical  approach  is  the  Gibbs  energy 
minimization  algorithm  (Greenberg  1986,  Harvie 
et  al.  1987,  Spencer  et  al.  1990).  Our  present  strat¬ 
egy  is  to  first  document,  by  using  the  FREZCHEM 
model,  where  convergence  problems  exists.  If  these 
problems  appear  to  be  intractable  using  our  present 
mathematical  algorithm,  we  plan  to  try  Gibbs 
energy  minimization. 
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PROGRAM  FREZCHEM 

C  FREZCHEM:  This  is  a  chemical  thermodynamic  model  for  aqueous 
C  solutions  over  the  temperature  range  from  +25C  to  -60C.  This 
C  model  uses  the  Pitzer  equations  for  activity  coefficients  and 
C  the  activity  of  water  and  is  valid  to  high  ionic  strengths 
C  (20  M) .  The  model  considers  the  precipitation-dissolution  of 
C  the  chloride  and  sulfate  salts  of  calcium,  magnesium,  sodium, 

C  and  potassium. 

REAL  IX, INIT(20) ,I,K 
DIMENSION  FINAL (20) 

CHARACTER* 50  TITLE 

DOUBLE  PRECISION  AKI, AKII , ISALGM, SALGM, DELGM, SALH20, X 
DOUBLE  PRECISION  M,M0,H20 

COMMON  /A/  MAXCAT,MAXAN,AKI (0:20) , AKII (0:20) ,BET0 (10,11:20) , 
xBETl  (10,11:20) ,BET2 (10, 11: 20) ,0(10,11:20) ,CPHI (10,11:20) , 
xPSI (20,20,20) , THETA (20, 20) ,A, Z(30) 

COMMON  /B/  X(60) ,GX(60) ,K(20:60) ,AH20,DEL(5) , 
xIX(30) , TITLE, ITER, AX(60) , TEMP, I, PHI 

OPEN  (UNIT  =  2,  FILE  =  "FrData") 

C  Read  in  or  initialize  run  parameters. 

PRINT*,  'TITLE  =' 

READ(*,1)  TITLE 
1  FORMAT  (A50) 

PRINT*,  'FREEZE (1)  OR  EVAPORATION (2)  SCENARIO  ?' 

READ*,  PATH 

PRINT*,  'SODIUM  (M/KG)  =' 

READ*,  IX  (1) 

PRINT*,  'POTASSIUM  (M/KG)  =' 

READ*,  IX  (2) 

PRINT*,  'CALCIUM  (M/KG)  =' 

READ*,  IX  (3) 

PRINT*,  'MAGNESIUM  (M/KG)  =' 

READ*,  IX  (4) 

PRINT*,  'CHLORIDE  (M/KG)  =' 

READ*,  IX (11) 

PRINT*,  'SULFATE  (M/KG)  =' 

READ*,  IX (12) 

PRINT*,  'INITIAL  TEMPERATURE  (K)  =' 

READ*,  TEMP 
IF  (PATH.EQ.l)  THEN 

PRINT*,  'FINAL  TEMPERATURE  (K)  =’ 

READ*,  FINTEM 

PRINT*,  'TEMPERATURE  DECREMENT  (K)  =' 

READ*,  DELTEM 

ELSE 

PRINT*,  'FINAL  WATER  (G)  =  ' 

READ*,  FINWAT 

PRINT*,  'WATER  DECREMENT  (G)  =  ' 

READ*,  DELWAT 

END  IF 

DATA  DEL  /  1.25,  1.10,  1.01,  1.001,  1.0001  / 

X(l)  =  IX(1) 

X(2)  =  IX(2) 

X(3)  =  IX(3) 

X(4)  =  IX(4) 

MAXCAT  =  4 
X(ll)  =  IX(ll) 
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X(12)  =  IX(12) 

MAXAN  =  12 
X(30)  =  1000 
X(31)  =  0 
ISALGM  =  0 
3  ITER  =  0 

CALL  PARAMETER 
5  CALL  PITZER 

C  Test  for  supersaturation  with  respect  to  ice  phase  and  freeze 
C  if  necessary. 

IF  (AH20  ,GT.  K{31))  THEN 
IF  (PATH.EQ.2)  THEN 

WRITE (2,*) 'THEORETICALLY  THIS  SOLUTION  SHOULD  FREEZE.' 
WRITE  (2,*)  'THEREFORE,  CAN  NOT  EVALUATE  THIS  USING  AN 
X  EVAPORATIVE  SCENARIO.' 

FINTEM=TEMP 
DELTEM=1 . 0 
PATH=1 

END  IF 
MO  =  0 
M  =  0 

H20  =  X(30)  +  X(31) 

DO  10  J  =  1,  25 

MO  =  MO  +  X(J) 

10  CONTINUE 

DO  70  KA  =  1,  5 

20  DO  30  J  =  1,  25 

X(J)  =  X(J)  *  DEL(KA) 

30  CONTINUE 

DO  40  J  =  32,  60 

X(J)  =  X(J)  *  DEL(KA) 

40  CONTINUE 

CALL  PITZER 

IF  (AH20  .GT.  K(31))  GOTO  20 
DO  50  J  =  1,  25 

X(J)  =  X(J)  /  DEL(KA) 

50  CONTINUE 

DO  60  J  =  32,  60 

X(J)  =  X(J)  /  DEL(KA) 

60  CONTINUE 

70  CONTINUE 

DO  80  J  =  1,  25 

M  =  M  +  X(J) 

80  CONTINUE 

X(30)  =  X(30)  *  MO  /  M 
X(31)  =  H20  -  X(30) 

END  IF 


C  Test  for  mineral  saturation  and  equilibrate  solutions. 


200 

210 


DO  210  J  =  1,  20 

INIT{J)  =  X(J) 


CONTINUE 
CALL  PITZER 
ITER  =  ITER  +  1 


IF 

((IX(3) 

.GT. 

0) 

.AND. 

(IX(12) 

.GT 

IF 

((IX(4) 

.GT. 

0) 

.AND. 

(IX(12) 

.GT 

IF 

((IX(1) 

.GT. 

0) 

.AND. 

(IX(12) 

.GT 

IF 

((IX(2) 

.GT. 

0) 

.AND. 

(IX  (12) 

.GT 

IF 

((IX(2) 

.GT. 

0) 

•  AND. 

(IX(4) 

.GT. 

0))  CALL  CAS04 
0 ) )  CALL  MGS04 
0))  CALL  NA2S04 
0))  CALL  ARCANITE 
0)  .AND.  (IX(12)  .GT.  0)) 
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X 

CALL 

PICROMERITE 

IF 

((IX(2) 

.GT.  0) 

.AND. 

(IX(ll) 

.GT.  0))  CALL 

SYLVITE 

IF 

((IX(1) 

.GT.  0) 

.AND. 

(IX  (ll) 

.GT.  0))  CALL 

NACL 

IF 

((IX(4) 

.GT.  0) 

.AND. 

(IXdl) 

.GT.  0))  CALL 

MGCL2 

IF 

((IX(3) 

.GT.  0) 

.AND. 

(IX (11) 

.GT.  0))  CALL 

ANTARCTICITE 

IF 

((IX(2) 

.GT.  0) 

.AND. 

(IX(4) 

.GT.  0)  .AND. 

(IX(ll)  .GT. 

0)) 

X 

CALL 

CARNALLITE 

IF 

((IX(3) 

.GT.  0) 

.AND. 

{IX(4) 

.GT.  0)  .AND. 

(IX(ll)  .GT. 

0)) 

X 

CALL 

TACHYHYDRITE 

IF 

(ITER  .: 

EQ.  300) 

THEN 

WRITE  (2,*) 

'THIS 

RUN  IS  NOT  CONVERGING  . 

AFTER  300 

XINTERATIONS ' 

GOTO  270 


END  IF 


C  Calculate  water  associated  with  hydrated  salts. 

SALH2O=2*X(32)+6*X(35)+6*X(36)+8*X(37)+12*X(38)+6*X{39)+12*X(40) 
SALH2O=SALH2O+10*X (41) +6*X (43) +7*X (44) +6*X (46) 

SALGM  =  SALH20  *  X(30)  •*  18.0153  /  1000 
DELGM  =  SALGM  -  I SALGM 
X(30)  =  X(30)  -  DELGM 
DO  220  J  =  1,  25 

X(J)  =  X(J)  *  (X(30)  +  DELGM)  /  X(30) 

220  CONTINUE 

DO  230  J  =  32,  60 

X(J)  =  X(J)  *  (X(30)  +  DELGM)  /  X(30) 

230  CONTINUE 

I SALGM  =  SALGM 

C  Test  for  mineral  equilibria. 

235  DO  240  J  =  1,  20 

FINAL (J)  =  X(J) 

IF  (FINAL (J)  .EQ.  0)  GOTO  240 

IF  ( (ABS (FINAL (J) -INIT(J) ) /FINAL (J) ) *1000  .GT.  1)  GOTO  200 
240  CONTINUE 

C  Test  for  ice-water  equilibrium 
CALL  PITZER 

IF  (AH20^K(31)  .GT.  0.00005)  GOTO  5 
C  Print  Results. 


260  CALL  PPRINT 

C  Update  Temperature  or  Water  Content 

IF  (PATH.EQ.l)  THEN 

TEMP  =  TEMP  -  DELTEM 
IF  (TEMP  .GE,  FINTEM)  GOTO  3 

ELSE 

X(30)=X(30)-DELWAT 
TH20  =  SALGM+X(30) 

IF  (TH20.lt. FINWAT)  GOTO  270 
DO  263  IJK=1,25 

X(IJK)=X(IJK) * (X(30)+DELWAT) /X(30) 
263  CONTINUE 

DO  265  IJK=32,60 

X(IJK)=X(IJK) * (X(30)+DELWAT) /X(30) 
265  CONTINUE 

GO  TO  3 
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END  IF 
270  END 

C - 

C - 


SUBROUTINE  PARAMETER 

C  This  subroutine  calculates  the  parameters  and  their  temperature 

C  dependence  for  the  Pitzer  equations  and  mineral  solubility  products. 

REAL  IX, I, K 

CHARACTER* 50  TITLE 

DOUBLE  PRECISION  D(64,6)',  P  (64)  ,AKI,AKII,T,X 

COMMON  /A/  MAXCAT,MAXAN,AKI (0:20) ,AKII (0:20) ,BET0 (10,11:20) , 
xBETl (10,11:20) ,BET2 (10, 11 : 20) , C (10, 11 : 20) , CPHI (10,11:20) , 
xPSI(20,20,20)  ,THETA(20,-20)  ,A,Z(30) 

COMMON  /B/  X(60) ,GX(60) ,K(20:60) ,AH20,DEL (5) , 
xIX(30) , TITLE, ITER, AX(60) , TEMP, I, PHI 

DATA  AKI  /  1.925154014814667,  -0.060076477753119, 

X  -0.029779077456514,  -0.007299499690937, 

X  0.000388260636404,  0.000636874599598, 

X  0.000036583601823,  -0.000045036975204, 

X  -0.000004537895710,  0.000002937706971, 

X  0.000000396566462,  -0.000000202099617, 

X  -0.000000025267769,  0.000000013522610, 

X  0.000000001229405,  -0.000000000821969, 

X  -0.000000000050847,  0.000000000046333, 

X  0.000000000001943,  -0.000000000002563, 

X  -0.000000000010991  / 

DATA  AKII  /  0.628023320520852,  0.462762985338493, 

X  0.150044637187895,  -0.028796057604906, 

X  -0.036552745910311,  -0.001668087945272, 

X  0.006519840398744,  0.001130378079086, 

X  -0.000887171310131,  -0.000242107641309, 

X  0.000087294451594,  0.000034682122751, 

X  -0.000004583768938,  -0.000003548684306, 

X  -0.000000250453880,  0.000000216991779, 

X  0.000000080779570,  0.000000004558555, 

X  -0.000000006944757,  -0.000000002849257, 

X  0.000000000237816  / 

DATA  Z  /  1,1, 2, 2, 1,0, 0,0, 0,0, -1,-2, -1,-1, -2, 

X  0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0,0  / 

C  Data  Stored  in  D(I,J)  are  from  Spencer  et  al.  (1990) 

DATA  ((D(I, J), J=1,6),I=1,64)  / 

X  86.6836498,  0.0848795942,  -8 . 8878515D-5, 

X  4.88096393D-8,  -1 . 32731477D3,  -17.6460172, 

X  7.87239712,  -8 . 3864096D-3,  1 . 44137774D-5, 

X  -8.7820301D-9,  -496.920671,  -0.820972560, 

X  866.915291,  0.606166931,  -4 . 8048921D-4, 

X  1.88503857D-7,  -1 . 70460145D4,  -167.171296, 

X  1.70761824,  2 . 32970177D-3,  -2. 46665619D-6, 

X  1.21543380D-9,  -1.35583596,  -0.387767714, 

X  26.5718766,  9 . 92715099D-3,  -3. 62323330D-6, 

X  -6.28427180D-11,  -755.707220,  -4.67300770, 
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X  1. 69742977D3,  1.22270943,  -9 . 99044490D-4, 

X  4.04786721D-7,  -3.2868442204,-328.813848, 

X  -3.27571680,  -1 . 27222054D-3,  4 . 71374283D-7, 

X  l.li62507D-ll,  90.7747666,  0.580513562, 

X  -5.62764702D1,  -3 . 00771997D-2,  1 . 05630400D-5, 

X  3.3331626D-9,  1.1173034903,  1 . 0666474301, 

X  3.4787,  -1.54170-2,  3.17910-5,  0,  0,  0, 

X  2.6423165501,  2.469229930-2,  -2.482985100-5, 

X  1.224218640-8,  -4 . 1809842702,  -5.35350322, 

X  3.1385291302,  2.617690990-1,  -2.462684600-4, 

X  1.157647870-7,  -5.5313338103,  -6.2161686201, 

X  -3.1843252504,  -2.8671035801,  2.788928380-2, 

X  -1.32797050-5,  5.2403295805,  6.4077039603, 

X  5.95320-2,  -2.499490-4,  2.418310-7,  0,  0,  0, 

X  -3.3248633003,  -2.9297353,  2 . 8024367O-3, 

X  -1.3168830-6,  5.5395852704,  6.6666036902, 

X  -3.5740616003,  -3.0011206,  2.736609500-3, 

X  -1.2191710-6,  6. 0971648204,  7.1161312002, 

X  3.6852047802,  3.162439950-1,  -2.953727600-4, 

X  1.354911040-7,  -6. 2260791303,  -7.3584409401, 

X  4.0790879701,  8.269066750-3,  0,0,  -1.4184299803,-6.74728848, 

X  -1.3166965101,  2.357932390-2,  0,  0,  2.0671259403,  0, 

X  -1.880-2,  0,  0,  0,  0,  0, 

X  1.500-1,  0,  0,  0,  0,  0, 

X  3.00,  0,  0,  0,  0,  0, 

X  -1.2939928702,  4 . 00431027O-1,  0,  0,  0,  0, 

X  5.4000784903,  4.90576884,  -4.804897500-3, 

X  2.311269940-6,  -8.8066414604,  -1.0883956503, 

X  2.78730869,  4.300774400-3,  0,  0,  0,  0, 

X  -5.8862365302,  -5.05522800-1,  4.82776570-4, 

X  -2.30298380-7,  1.0200201604,  1.1730380802, 

X  -18.2266741,  -3.690384700-3,  0, 

X  0,  612.415011,  3.02994981, 

X  6.48108127,  1.468034680-3,  0, 

X  0,  -204.354019,  -1.09448043, 

X  3.481200-2,  0,  0,  0,  -8.21660,  0, 

X  5.00-2,  0,  0,  0,  0,  0, 

X  -7.63980,  -1.29900-2,  1.10600-5,  0,  0,  1.8475, 

X  -1.200-2,  0,  0,  0,  0,  0, 

X  7.00-2,  0,  0,  0,  0,  0, 

X  -3.1098700-2,  5.44647800-5,  0,  0,  1.99404210,  0, 

X  1.17505200-1,  0,  0,  0,  -4.19862001,  0, 

X  2.365710,  -4.5400-3,  0,  0,  -2.8494002,  0, 

X  -5.9300-2,  2.542800-4,  0,  0,  -1.3439001,  0, 

X  1.16700-1,  0,  0,  0,  0,  0, 

X  5.03622300-2,  -8.7508200-6,  0,  0,  -2.89909001,  0, 

X  -1.3679157,  4.240166530-3,  0,  0,  0,  0, 

X  5.31274136,  -6.34242480-3,  0,  0,  -9.8311384702,  0, 

X  4.1579022001,  1.303773120-2,  0,0,  -9.8165852602,  -7.4061986, 

X  7.000-2,  0,  0,  0,  0,  0, 

X  4.0209775,  1.12860050-3,  0,  0,  -1.0116926002,  -7.0607980-1, 

X  -2.1248150-1,  2.84698330-4,  0,  0,  3.7561961401,  0, 

X  -1.800-2,  0,  0,  0,  0,  0, 

X  -1.8391580-1,  1.4294440-4,  0,  0,  3.26301,  0, 

X  -1.6291734103,  -1.51940390,  1.452496790-3, 

X  -6.94275050-7,  2.2601274304,  333.075506, 

X  7.87506039303,  11.69118490,  -1.71837890-2, 

X  1.243955430-5,  -9.331479004,  -1.728746103, 

X  -1.222255104,  -9.8806459,  8.466850830-3, 

X  -3.44591170-6,  2.0982396505,  2.4232852803, 

X  -1.4747774501,  0,  0,  0,  3.2640949603,  0, 

X  3.6528311502,  -7.15782570-1,  0,  0,  -4.4875339104,  0, 

X  9.1483900103,  8.22348745,  -8.12887590-3, 
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730 


3.95552403D-6,  -1 . 54040868D5,  -1 . 83624247D3, 
1.42290062D5,  1 . 61973105D2,  -1 . 95332071D-1, 
1.17636119D-4,  -2 . 04059847D6,  -2. 97464810D4, 
7.52225099D2,  1.175846530-1,0,0,-2.4322390904,- 
2.2780197603,  6.493616160-1,0,0,-6.2307512304,- 
2.5500889605,  2.4453224002,  -2.488078760-1, 
1.224252360-4,  -4.0298834206,  -5. 1866860404, 
-4.4570217101,  2.320237900-1,  -7.149356920-4, 
5.326582150-7,  -4.2481792303,  8.59110245, 
8.0377791801,  -1.3880690-1,  0,  0,  0,  0, 

8. 4472805004,  7.6844338701,  -7.48258160-2, 
3.518060850-5,  -1.388185206,  -1,702677804, 
-3.963563203,  -5.8114490,  7.597994620-3, 
-4.65717370-6,  3.9345489304,  8.7959842302, 
2.3916646101,  -9.36807440-2,  0,  0,  0,  0, 
2.02443627,  -1.80515860-2,  0,  0,  0,  0, 
-3.5334625104,  -2.7680999101,  2.398560210-2, 
-1.02018450-5,  6.3756198805,  6.9473778703, 
-1.7520350901,  2.7937130-2,  0,  0,  0,  0  / 


T=TEMP 

00  730  J  =  1,  64 

P ( J) =0 ( J, 1) +0 ( J, 2) *T+0 (J, 3) *T**2+0 ( J,  4) *T**3+0 ( J, 
/T+0(J, 6)*LOG(T) 

CONTINUE 


A  =  P(l) 
BET0(1,11) 

_ 

P(2) 

BETl (1,11) 

= 

P(3) 

CPHI (1,11) 

3S 

P(4) 

BETO (2,11) 

5= 

P(5) 

BETl (2, 11) 

= 

P(6) 

CPHI (2, 11) 

s- 

P(7) 

BETO  (3, 11) 

5= 

P(8) 

BETl (3, 11) 

= 

P(9) 

CPHI (3, 11) 

= 

P(10) 

BETO (4,11) 

P(ll) 

BETl (4, 11) 

=: 

P(12) 

CPHI  (4,11) 

=: 

P(13) 

BETO (1,12) 

P(14) 

BETl (1,12) 

= 

P(15) 

CPHI (1,12) 

= 

P(16) 

BETO (2,12) 

= 

P(17) 

BETl (2,12) 

= 

P(18) 

CPHI (2, 12) 

= 

P(19) 

BETO  (3,12) 

= 

P(20) 

BETl (3,12) 

= 

P(21) 

BET2 (3, 12) 

= 

P(22) 

BETO (4, 12) 

= 

P(23) 

BETl (4,12) 

= 

P(24) 

CPHI (4,12) 

= 

P(25) 

THETA (1,2) 

= 

P(26) 

PSI  (1,2,11) 

=  P(27) 

PSI  (1,2,12) 

=  P(28) 

THETA (1,  3) 

= 

P(29) 

PSI  (1,3,11) 

=  P(30) 

PSI(1,3,12) 

=  P(31) 

THETAd,  4) 

= 

P(32) 

PSI (1,4, 11) 

=  P(33) 

PSI (1,4, 12) 

=  P(34) 

THETA (2, 3) 

_ 

P(35) 

PSI  (2, 3, 11) 

=  P(36) 

THETA (2, 4) 

= 

P(37) 

PSI (2, 4, 11) 

=  P(38) 

.2199007602, 

.9543889102, 
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PSI(2,4,12)  =  P(39) 
THETA(3,4)  =  P(40) 
PSI(3,4,11)  =  P(41) 
PSI(3,4,12)  =  0.024 
THETA(11,12)  =  P(42) 
PSI(11,12,1)  =  P(43) 
PSI(11,12,2)  =  P(44) 
PSI(11,12,3)  =  P(45) 
PSI(11,12,4)  =  P(46) 
K(34)  =  EXP(P(47) ) 
K(31)  =  EXP(P(48)  ) 
K(32)  =  EXP(P(49)  ) 
K{22)  =  EXP(P(50)  ) 
K(23)  =  EXP (P  (51) ) 
K{33)  =  EXP (P  (52) ) 
K(35)  =  EXP(P(53)  ) 
K(36)  =  EXP(P(54)  ) 
K(37)  =  EXP(P(55)  ) 
K(38)  =  EXP(P(56)  ) 
K(39)  =  EXP (P (57) ) 
K(40)  =  EXP(P(58)  ) 
K(41)  =  EXP(P(59)  ) 
K(42)  =  EXP(P(60) ) 
K(43)  =  EXP(P(61)  ) 
K(44)  =  EXP(P(62)  ) 
K(45)  =  EXP(P(63)  ) 
K(46)  =  EXP(P(64)  ) 
RETURN 
END 


C - 

C - 

SUBROUTINE  PITZER 

C  This  subroutine  calculates  activity  coefficients  and  the  activity 
C  of  water  using  the  Pitzer  equations. 

REAL  I,IX,K 

CHARACTER* 50  TITLE 

DOUBLE  PRECISION  AKI,AKII,X 

DIMENSION  SUMCA (10) ,  SUMCAT(IO),  SUMAN(IO),  SUMZ (10) , 

X  SUMAC (11: 20) ,  SUMAA(11:20) ,  SUMCC (11 : 20) ,  SUMK(11:20), 

X  BPHKIO, 11:20),  B(10, 11:20),  BPRIME  (10, 11 : 20) , 

X  PHIPHI(20,20),  PHIPRI (20,20),  PHIIJ(20,20) 

COMMON  /A/  MAXCAT,MAXAN,AKI(0:20),AKII(0:20),BETO(10, 11:20), 
xBETl(10,ll:20),BET2(10,ll:20),C(10,ll:20),CPHI(10,ll:20), 

XPSI (20,20,20) , THETA (20, 20) ,A,Z(30) 

COMMON  /B/  X(60),GX(60),K(20:60),AH2O,DEL(5), 

XIX(30) , TITLE, ITER,AX(60) , TEMP, I, PHI 

ZZ  =  0 

SM  =  0 

1  =  0 

DO  280  J  =  1,  20 

ZZ  =  ZZ  +  X(J)  *  ABS(Z(J)) 

1  =  1+  .5  *  (X(J)  *  Z(J)**2) 

SM  =  SM  +  X(J) 

280  CONTINUE 

DATA  SUMCA  /  10*0  / 
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DATA  SUMCAT  /  10*0  / 

DATA  SUMAN  /  10*0  / 

DATA  SUMZ  /  10*0  / 

DATA  SUMAC  /  10*0  / 

DATA  SUMAA  /  10*0  / 

DATA  SUMCC  /  10*0  / 

DATA  SUMK  /  10*0  / 

ALPHA  =  2  *  SORT (I) 

ALPHAl  =  1.4  *  SORT (I) 

ALPHA2  =  12  *  SQRT(I) 

G1  =  2* (1- (1+ALPHAl) *EXP (-ALPHA1) ) /ALPHAl **2 

G2  =  2* (1- (1+ALPHA2) *EXP (-ALPHA2) ) /ALPHA2**2 

GPRIl  =  -2* (1- (l+ALPHAl+ALPHAl**2/2) *EXP (-ALPHAl) ) /ALPHA1**2 

GPRI2  =  -2*(1-(1+ALPHA2+ALPHA2**2/2)*EXP(-ALPHA2) )/ALPHA2**2 

G  =  2* (1- (1+ALPHA) *EXP (-ALPHA) ) /ALPHA**2 

GPRIME  =  -2* (1- (l+ALPHA+ALPHA**2/2) *EXP (-ALPHA) ) /ALPHA**2 
DO  300  J  =  1,  MAXCAT 

DO  290  IK  =  11,  MAXAN 

IF  (Z(J)*ABS(Z(IK))  .EQ.  4)  THEN 

BPHI ( J, IK) =BET0 ( J, IK) +BET1 ( J, IK) *EXP (-ALPHAl) + 
X  BET2 (J, IK)*EXP(-ALPHA2) 

B ( J, IK) =BET0 ( J, IK) +BET1 ( J, IK) *G1+BET2 ( J, IK) *G2 
BPRIME ( J, IK) =BET1 ( J,  IK) *GPRI1/I+BET2 ( J,  IK) * 

X  GPRI2/I 

ELSE 

BPHI ( J, IK) =BET0 ( J, IK) +BET1 ( J, IK) *EXP (-ALPHA) 

B ( J, IK) =BET0 ( J, IK) +BET1 ( J, IK) *G 
BPRIME ( J, IK) =BET1 ( J, IK) *GPRIME/I 

END  IF 

C ( J, IK) =CPHI ( J, IK) / (2*SQRT (Z ( J) *ABS (Z (IK) ) ) ) 

290  CONTINUE 

300  CONTINUE 

DO  320  J  =  1,  MAXCAT  -  1 

DO  310  IK  =  J  +  1,  MAXCAT 

IF  ((Z(J)  .EQ.  0)  .OR.  (Z(IK)  ,EQ.  0))  GOTO  310 
CALL  INTERACTION (Z ( J) , Z (IK) , I, A, PHIPHI (J, IK) , 

X  PHIPRI ( J,  IK) ,PHIIJ( J,  IK) , THETA ( J, IK) , AKI, AKII) 

310  CONTINUE 

320  CONTINUE 

DO  340  J  =  11,  MAXAN  -  1 

DO  330  IK  =  J  +  1,  MAXAN 

IF  ((Z(J)  .EQ.  0)  .OR.  (Z(IK)  .EQ.  0))  GOTO  330 
CALL  INTERACTION (Z ( J) , Z (IK) , I, A, PHIPHI ( J, IK) , 

X  PHIPRI (J,IK) ,PHIIJ(J, IK) ,THETA(J,IK) , AKI, AKII) 

330  CONTINUE 

340  CONTINUE 

C  Calculation  of  summation  terms  for  F  and  PHI. 

SCATON  =  0 
SUBSUM  =  0 
SANON  =  0 
SUMCAF  =  0 
SUMANF  =  0 

DO  370  J  =  1,  MAXCAT  -  1 

DO  360  IK  =  J  +  1,  MAXCAT 
DO  350  L  =  11,  MAXAN 

SUBSUM  =  SUBSUM  +  PSI(J,IK,L)  *  X(L) 

350  CONTINUE 

SCATON=SCATON+ (SUBSUM+PHIPHI ( J, IK) ) *X ( J) *X (IK) 
SUMCAF=SUMCAF+PHIPRI ( J, IK) *X ( J) *X (IK) 
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SUBSUM  =  0 

360  CONTINUE 

370  CONTINUE 

SUBSUM  =  0 

DO  400  J  =  11,  MAXAN  -  1 

DO  390  IK  =  J  +  1,  MAXAN 

DO  380  L  =  1,  MAXCAT 

SUBSUM=SUBSUM+PSI ( J, IK, L) *X (L) 

380  CONTINUE 

SANON=SANON+ (SUBSUM+PHIPHI ( J, IK) ) *X ( J) *X (IK) 
SUMANF=SUMANF+PHIPRI ( J, IK) *X ( J) *X (IK) 

SUBSUM  =  0 

390  CONTINUE 

400  CONTINUE 
SUMB  =  0 
SUMPHI  =  0 

DO  420  J  =  1,  MAXCAT 

DO  410  IK  =  11,  MAXAN 

SUMB=SUMB+X ( J) *X (IK) *BPRIME ( J, IK) 

SUMPHI=SUMPHI+X { J) *X (IK) * (BPHI ( J, IK) +ZZ*C ( J, IK) ) 
410  CONTINUE 

420  CONTINUE 

F=-A* (SORT (I) / (1+1.2*SQRT(I) )+2*LOG(l+1.2*SQRT(I) ) /1.2)+SUMB+ 
X  SUMCAF+SUMANF 

PHI=1+ (2/SM) * ( (-A*I**1.5) / (1+1.2*SQRT(I) ) +SUMPHI+SCATON+ 

X  SANON) 

AH20  =  EXP (-PHI*SM/55. 50837) 

C  Calculation  of  terms  for  activity  coefficients  (gamma) . 

SUM  =  0 

DO  450  J  =  1,  MAXCAT  -  1 

DO  440  IK  =  J  +  1,  MAXCAT 
DO  430  L  =  11,  MAXAN 

PSI(IK,J,L)  =  PSI(J,IK,L) 

PHIIJ(IK,J)  =  PHIIJ(J, IK) 

430  CONTINUE 

440  CONTINUE 

450  CONTINUE 

DO  480  L  =  11,  MAXAN  -  1 

DO  470  IK  =  L  +  1,  MAXAN 

DO  460  J  =  1,  MAXCAT 

PSI(IK,L,J)  =  PSI(L,  IK,  J) 

PHIIJ(IK,L)  =  PHIIJ(L,IK) 

460  CONTINUE 

470  CONTINUE 

480  CONTINUE 

DO  500  J  =  1,  MAXCAT 

DO  490  IK  =  11,  MAXAN 

SUMCA(J)  =  SUMCA(J)+X(IK)*(2*(B(J,IK))+ZZ*C(J,IK)) 
490  CONTINUE 

500  CONTINUE 

DO  530  J  =  1,  MAXCAT 

DO  520  IK  =  1,  MAXCAT 

IF  (J  .EQ.  IK)  GOTO  520 
DO  510  L  =  11,  MAXAN 

SUM=SUM+X (L) *PSI ( J, IK, L) 

510  CONTINUE 

SUMCAT ( J) =SUMCAT ( J) +X (IK) * (SUM+2*PHIIJ( J, IK) ) 

SUM  =  0 

520  CONTINUE 

530  CONTINUE 

DO  560  J  =  1,  MAXCAT 
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DO  550  IK  =  11,  MAXAN-1 

DO  540  L  =  IK  +  1,  MAXAN 

SUMAN(J)  =  SUMAN(J)+X{IK)*X(L)*PSI(IK,L, J) 

540  CONTINUE 

550  CONTINUE 

560  CONTINUE 
SUM  =  0 

DO  580  J  =  1,  MAXCAT 

DO  570  IK  =  11,  MAXAN 

SUM  =  SUM  +  X(J)*X(IK)*C(J, IK) 

570  CONTINUE 

580  CONTINUE 

DO  590  J  =  1,  MAXCAT 

SUMZ(J)  =  SUMi^ABS(Z(J)  ) 

590  CONTINUE 

DO  600  J  =  1,  MAXCAT 

GX ( J) =EXP (Z ( J) **2*F+SUMCA ( J) +SUMCAT ( J) +SUMAN ( J) +SUMZ ( J) ) 

600  CONTINUE 
SUM  =  0 

DO  620  IK  =  11,  MAXAN 

DO  610  J  =  1,  MAXCAT 

SUMAC (IK)  =  SUMAC{IK)+X(J)*(2*B(J,IK)+ZZ*C(J,IK)) 

610  CONTINUE 

620  CONTINUE 

DO  650  IK  =  11,  MAXAN 

DO  640  L  =  11,  MAXAN 

IF  (IK  .EQ.  L)  GOTO  640 
DO  630  J  =  1,  MAXCAT 

SUM  =  SUM+X(J)*PSI(IK,L, J) 

630  CONTINUE 

SUMAA (IK) =SUMAA (IK) +X (L) * (2*PHII J (IK, L) +SUM) 

SUM  =  0 

640  CONTINUE 

650  CONTINUE 

DO  680  IK  =  11,  MAXAN 

DO  670  J  =  1,  MAXCAT  -  1 

DO  660  L  =  J  +  1,  MAXCAT 

SUMCC (IK) =SUMCC (IK) +X ( J) *X (L) *PSI ( J, L, IK) 

660  CONTINUE 

670  CONTINUE 

680  CONTINUE 
SUM  =  0 

DO  700  J  =  1,  MAXCAT 

DO  690  IK  =  11,  MAXAN 

SUM  =  SUM+X(J)*X(IK)*C(J,IK) 

690  CONTINUE 

700  CONTINUE 

DO  710  IK  =  11,  MAXAN 

SUMK(IK)  =  SUM  *  ABS(Z(IK)) 

710  CONTINUE 

DO  720  IK  =  11,  MAXAN 

GX (IK) =EXP (Z (IK) **2*F+SUMAC (IK) +SUMAA (IK) +SUMCC (IK) +SUMK (IK) ) 
720  CONTINUE 
RETURN 
END 

C - 

C - 


SUBROUTINE  INTERACTION (Zl, Z2, I, A, PHIPHI, PHIPRI, PHIIJ, THETA, 
xAKI,AKII) 
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C  This  subroutine  calculates  the  higher-order  electrostatic 
C  interaction  terms  for  the  Pitzer  equations. 

REAL  I,  JX(2,2),  JXPRIME(2,2) 

DIMENSION  B (0:22) ,  D (0 : 22) , ZA (2) ,  XA(2,2) 

DOUBLE  PRECISION  AKI (0 : 20) , AKII (0 : 20) 

B(21)  =  0 
B(22)  =  0 
D(21)  =  0 
D(22)  =  0 
ZA(1)  =  Z1 
ZA(2)  =  Z2 

DO  770  JJ  =  1,  2 

DO  760  IJ  =  1,  2  ' 

XA ( JJ, I J) =6*ZA (IJ) *ZA( JJ) *A*SQRT (I) 

X=XA(JJ, IJ) 

IF  (X  .LT.  1)  THEN 

ZZ  =  ■4*X**0. 20-2.0 
DZ  =  0.80*X** (-0.80) 

DO  740  K  =  20,  0,  -1 

B(K)  =  ZZ*B(K+l)-B(K+2)+AKI(K) 

D(K)  =  B(K+l)+ZZ*D(K+l)-D(K+2) 

740  CONTINUE 

ELSE 

ZZ  =  (40/9)*X**(-0.10)-22/9 
DZ  =  (-40/90)*X** (-1.10) 

DO  750  K  =  20,  0,  -1 

B(K)  =  ZZ*B(K+l)-B(K+2)+AKII(K) 

D(K)  =  B(K+l)+ZZ*D(K+l)-D(K+2) 

750  CONTINUE 

END  IF 

JX(JJ,IJ)  =  ,25*X-1+.5*(B(0)-B(2)) 

JXPRIME(JJ,IJ)  =  .25+.5*DZ*(D(0)-D(2)) 

760  CONTINUE 

770  CONTINUE 

ETHETA=(Z1*Z2/ (4*1 ) ) * ( JX  (1, 2) - . 5* JX (1, 1 ) - . 5* JX (2, 2) ) 

ETHPRI=-ETHETA/I+ (Z1*Z2/ (8*I**2) ) * (XA(1, 2) *JXPRIME (1, 2) -. 5* 

X  XA(1,1)*JXPRIME(1,1)-.5*XA(2,2)*JXPRIME(2,2) ) 

PHIPHI=THETA+ETHETA+I*ETHPRI 
PHIIJ=THETA+ETHETA 
PHIPRI=ETHPRI 
RETURN 
END 

C - 

C - 

SUBROUTINE  PPRINT 
C  Subroutine  prints  the  results. 

REAL  MOLE,  K,  IX,  I 

CHARACTER*20  SPECIS(30),  SOLID (31: 60) 

CHARACTER* 50  TITLE 
DIMENSION  AX (60) 

DOUBLE  PRECISION  X 

COMMON  /B/  X(60),GX(60),K(20:60),AH20,DEL(5), 
xIX (30)  ,  TITLE, ITER, AX (60) , TEMP, I, PHI 
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DATA  SPECIS  /  'NA', 'K', 'CA', 'MG', 'H', 0,0, 0,0,0, 

X  'CL',  ’S04', 'OH', 'HC03', 'C03' , 0, 0, 0, 0, 0, 

X  'C02(AQ) ', 'CAS04', 'MGS04', 0,0, 0,0,0, 'C02 (ATM) ', 'H20(L) '  / 

DATA  SOLID  /  ' ICE ' , 'NACL. 2H20' , 'NACL* , ' KCL ' , ' CACL2 . 6H20' , 

X  'MGCL2.6H20', 'MGCL2.8H20', 'MGCL2 . 12H20' , 'KMGCL3. 6H20' , 

X  'CACL2.2MGCL2.12H20', 'NA2SO4.10H2O' , 'NA2S04', 'MGS04 . 6H20' , 

X  'MGS04.7H20', 'K2S04', 'MGS04 .K2S04 . 6H20' , 

X  0,0, 0,0, 0,0, 0,0, 0,0, 0,0, 0,0  / 

GX(22)  =  1 
GX(23)  =  1 
CALL  PITZER 
DO  780  J  =  1,  30 

AX(J)  =  X(J)  *  GX(J) 

780  CONTINUE 

WRITE (2, 785)  TITLE 
785  FORMAT (A50) 

WRITE  (2,*) 

WRITE  (2,787)  "Temp(K) ", "Ion. Str . "AH20", "Phi", "H20 (g) ", "Ice(g) " 
787  FORMAT  (IX,  A7,  A15,  AlO,  A13,  A14,  A13) 

WRITE  (2,790)  TEMP, I, AH20, PHI,X (30) , X (31 ) 

790  FORMAT  (IX,  G10.5,  IX,  5G13.5) 

WRITE  (2,*) 

WRITE  (2,795)  "Solution", "Initial", "Final" 

795  FORMAT  (IX,  A8,  A13,  All) 

WRITE  (2,797)  "SPECIES", "Conc. ", "Cone, ", "Act .Coef . ", "Activity", 

X  "Moles" 

797  FORMAT  (IX,  A7,  A12,  A13,  A17,  A12,  AlO) 

DO  810  J  =  1,  23 

MOLE  =  X(J)*X(30)/1000 
IF  (MOLE  .EQ.  0)  GOTO  810 

WRITE (2, 800)  SPECIS (J),  IX (J),  X(J),  GX(J),  AX(J),  MOLE 
800  FORMATdX,  All,  5G13.5) 

810  CONTINUE 

MOLE  =  X(30)  /  18.0153 
WRITE  (2,812)  SPECIS (30),  AH20,  MOLE 
812  FORMAT  (IX,  All,  42X,  G10.5,G12.5) 

WRITE  (2,*) 

WRITE  (2,815)  "Solid", "Equil." 

815  FORMAT  (IX,  A5,30X,A7) 

WRITE  (2,817)  "SPECIES", "Moles", "Constant" 

817  FORMAT  (IX,  A7,  A21,  A16) 

DO  840  J  =  31,  46 

MOLE  =  X(J)*X(30)/1000 

IF  (J  .EQ,  31)  MOLE  =  X(31)  /  18.0153 

WRITE  (2,830)  SOLID (J),  MOLE,  K(J) 

830  FORMAT  (IX,  A20,2G13.5) 

840  CONTINUE 

WRITE  (2,*) 

WRITE  (2,845)  "CaS04  (ionpair)",  K(22) 

845  FORMAT  (IX,  A15,  G31.5) 

WRITE  (2,850)  "MgS04  (ionpair)",  K(23) 

850  FORMAT  (IX,  A15,  G30.5) 

WRITE  (2,*) 

WRITE  (2,*)  "  Iterations  =",  ITER 
WRITE  (2,*) 

WRITE  (2,*) 

RETURN 

END 

C - 

C - 
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SUBROUTINE  NACL 

C  Subroutine  for  NaCl  equilibrium. 

REAL  IX, I, K 
CHARACTER*50  TITLE 

DOUBLE  PRECISION  NAT, CLT, KPRIME, X, DELTA 

COMMON  /B/  X(60) ,GX(60) ,K(20:60) ,AH20,DEL(5) , 
xIX(30) , TITLE, ITER, AX (60) , TEMP, I, PHI 

DELTA=0 

HHALIT  =  K(32)  /  AH20**2 
HALITE  =  K(33) 

NAT  =  X(l)  +  X(32)  +  X(33) 

CLT  =  X(ll)  +  X(32)  +  X(33) 

IF  (HHALIT  .GT.  HALITE)  KPRIME  =  K (33 ) / (GX (1 ) *GX (11 ) ) 

IF  (HALITE  .GT.  HHALIT)  KPRIME  =  K (32) / (GX (1 ) *GX (11 ) *AH20**2 ) 
IF  (NAT*CLT-KPRIME  .GT.  0)  GOTO  860 
855  X(l)  =  NAT 
X(ll)  =  CLT 
X(32)  =  0 
X(33)  =  0 
GOTO  880 

860  DECREM  =  MIN (NAT, CLT) 

DO  870  KA  =  2,  5 

865  DELTA  =  DELTA+DECREM* (DEL (KA)-l) 

IF  (DELTA  .LT.  0)  GOTO  855 
X(l)  =  NAT  -  DELTA 
X(ll)  =  CLT  -  DELTA 
IF  (KA  .EQ.  5)  CALL  PITZER 
HHALIT  =  K(32)  /  AH20**2 
HALITE  =  K(33) 

IF  (HHALIT  .GT.  HALITE)  KPRIME  =  K (33) / (GX (1) *GX (11) ) 

IF  (HALITE  .GT.  HHALIT)  KPRIME  =  K ( 32 ) / (GX (1 ) *GX (11 ) * 

X  AH20**2) 

IF  (X(1)*X (11) -KPRIME  .GT.  0)  GOTO  865 
DELTA  =  DELTA-DECREM* (DEL(KA)-l) 

870  CONTINUE 

X(l)  =  NAT  -  DELTA 
X(ll)  =  CLT  -  DELTA 
IF  (HHALIT  .GT.  HALITE)  THEN 
X(33)  =  DELTA 
X(32)  =  0 

ELSE 

X(32)  =  DELTA 
X(33)  =  0 

END  IF 
880  RETURN 
END 


SUBROUTINE  SYLVITE 

C  Subroutine  for  KCl  equilibrium. 

REAL  IX, I, K 
CHARACTER*50  TITLE 

DOUBLE  PRECISION  KT, CLT, KPRIME, X, DELTA 
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COMMON  /B/  X(60),GX(60),K{20:60),AH20,DEL(5), 
xIX (30) ,  TITLE, ITER, AX (60) , TEMP, I, PHI 

DELTA=0 

KT  =  X(2)  +  X(34) 

CLT  =  X(ll)  +  X(34) 

KPRIME  =  K(34)  /  (GX (2) *GX (11) ) 

IF  (KT*CLT-KPRIME  .GT.  0)  GOTO  890 
885  X(2)  =  KT 

X(ll)  =  CLT 
X(34)  =  0 
GOTO  910 

890  DECREM  =  MIN (KT, CLT) 

DO  900  KA  =  2,  5 

895  DELTA  =  DELTA+DECREM* (DEL (KA) -1) 

IF  (DELTA  .LT.  0)  GOTO  885 
X(2)  =  KT  -  DELTA 
X(ll)  =  CLT  -  DELTA 
IF  (KA  .EQ.  5)  CALL  PITZER 
KPRIME  =  K(34)/(GX(2)*GX(11)) 

IF  (X (2) *X (11) -KPRIME  .GT.  0)  GOTO  895 
DELTA  =  DELTA-DECREM* (DEL(KA)-l) 

900  CONTINUE 

X(2)  =  KT-DELTA 
X(ll)  =  CLT-DELTA 
X(34)  =  DELTA 
910  RETURN 
END 

C - 

C - 


SUBROUTINE  ANTARCTICITE 

C  Subroutine  for  CaC12  equilibrium. 

REAL  IX, I, K 
CHARACTER*50  TITLE 

DOUBLE  PRECISION  CAT, CLT, KPRIME, X, DELTA 

COMMON  /B/  X(60) ,GX(60) ,K(20:60) ,AH20,DEL(5) , 
xIX(30) , TITLE, ITER,AX(60) , TEMP, I, PHI 

DELTA=0 

CAT  =  X(3)  +  X(35) 

CLT  =  X(ll)  +  2*X(35) 

KPRIME  =  K(35)/(GX(3)*GX(11)**2*AH20**6) 

IF  (CAT*CLT**2-KPRIME  .GT.  0)  GOTO  920 
915  X(3)  =  CAT 

X(ll)  =  CLT 
X(35)  =  0 
GOTO  940 

920  DECREM  =  MIN (CAT, CLT) 

DO  930  KA  =  2,  5 

925  DELTA  =  DELTA+DECREM* (DEL (KA) -1) 

IF  (DELTA  .LT.  0)  GOTO  915 

X(3)  =  CAT  -  DELTA 

X(ll)  =  CLT  -  2*DELTA 

IF  (KA  .EQ.  5)  CALL  PITZER 

KPRIME  =  K(35)/(GX(3)*GX(11)**2*AH20**6) 

IF  (X(3) *X(11) **2-KPRIME  .GT.  0)  GOTO  925 

DELTA  =  DELTA  -  DECREM* (DEL (KA) -1 ) 
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930  CONTINUE 

X(3)  =  CAT-DELTA 
X(ll)  =  CLT  -  2*DELTA 
X(35)  =  DELTA 
940  RETURN 
END 

C - 

C - 

SUBROUTINE  MGCL2 

C  Subroutine  for  MgC12  equilibrium. 

REAL  IX, I, K 
CHARACTER* 50  TITLE 

DOUBLE  PRECISION  MGT, CLT, KPRIME, X, DELTA, KPRIMl, KPRIM2, KPRIM3 

COMMON  /B/  X{60)  ,GX(60)-,K(20:60)  ,AH20,DEL(5)  , 
xIX ( 30)  ,  TITLE, ITER, AX (60) , TEMP, I, PHI 

DELTA=0 

MGT  =  X(4)+X(36)+X(37)+X(38) 

CLT  =  X(11)+2*X(36)+2*X(37)+2*X(38) 

KPRIMl  =  K(36)/(GX(4)*GX(11)**2*AH20**6) 

KPRIM2  =  K(37)/(GX(4)*GX(11)**2*AH20**8) 

KPRIM3  =  K(38)/(GX(4)*GX(11)**2*AH20**12) 

KPRIME=MIN (KPRIMl , KPRIM2 , KPRIM3) 

IF  (MGT*CLT**2-KPRIME  .GT.  0)  GOTO  960 
950  X(4)  =  MGT 

■'X(ll)  =  CLT 
X(36)  =  0 
X(37)  =  0 
X(38)  =  0 
GOTO  990 

960  DECREM  =  MIN (MGT, CLT) 

DO  980  KA  =  2,  5 

970  DELTA  =  DELTA+DECREM* (DEL (KA) -1) 

IF  (DELTA  .LT.  0)  GOTO  950 

X(4)  =  MGT  -  DELTA 

X(ll)  =  CLT-2*DELTA 

IF  (KA  .EQ.  5)  CALL  PITZER 

KPRIMl  =  K(36)/(GX(4)*GX(11)**2*AH20**6) 

KPRIM2  =  K(37)/(GX(4)*GX(11)**2*AH20**8) 

KPRIM3  =  K(38)/(GX(4)*GX(11)**2*AH20**12) 

KPRIME=MIN (KPRIMl, KPRIM2, KPRIM3) 

IF  (X(4)*X(11)**2-KPRIME  .GT.  0)  GOTO  970 
DELTA  =  DELTA  -  DECREM* (DEL (KA) -1 ) 

980  CONTINUE 

X(4)  =  MGT  -  DELTA 
X(ll)  =  CLT  -  2*DELTA 
IF  (KPRIMl  .EQ.  KPRIME)  THEN 
X(36)  =  DELTA 
X(37)  =  0 
X(38)  -  0 

END  IF 

IF  (KPRIM2  .EQ.  KPRIME)  THEN 
X(37)  =  DELTA 
X(36)  =  0 
X(38)  =  0 

END  IF 

IF  (KPRIM3  .EQ.  KPRIME)  THEN 
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990 


X(38)  =  DELTA 
X(36)  =  0 
X(37)  =  0 

END  IF 
RETURN 
END 


C 

C 


SUBROUTINE  CARNALLITE 

C  Subroutine  for  KMgC13.6H20  equilibrium. 

REAL  IX, I, K 
CHARACTER* 50  TITLE 

DOUBLE  PRECISION  KT,MGT, CLT, KPRIME, X, DELTA 

COMMON  /B/  X(60) ,GX(60)>K(20:60) ,AH20,DEL(5) , 
xIX (30)  ,  TITLE, ITER, AX (60) , TEMP, I, PHI 

DELTA=0 

KT  =  X(2)  +  X(39) 

MGT  =  X(4)  +  X(39) 

CLT  =  X(ll)  +  3*X(39) 

KPRIME  =  K(39)/(GX(2)*GX(4)*GX(11)**3*AH20**6) 

IF  (KT*MGT*CLT**3-KPRIME  .GT.  0)  GOTO  1005 
1000  X(2)  =  KT 

X(ll)  =  CLT 
X(4)  =  MGT 
X(39)  =  0 
GOTO  1030 

1005  DECREM  =  MIN (KT, MGT, CLT) 

DO  1020  KA  =  2,  5 

1010  DELTA  =  DELTA+DECREM* (DEL(KA)-l) 

IF  (DELTA  .LT.  0)  GOTO  1000 
X(2)  =  KT  -  DELTA 
X(4)  =  MGT  -  DELTA 
X(ll)  =  CLT  -  3*DELTA 
IF  (KA  ,EQ.  5)  CALL  PITZER 

KPRIME  =  K(39)/(GX(2)*GX(4)*GX(11)**3*AH20**6) 
IF  (X(2)*X(4)*X(11)**3-KPRIME  .GT.  0)  GOTO  1010 
DELTA  =  DELTA  -  DECREM* (DEL (KA) -1 ) 

1020  CONTINUE 

X(2)  =  KT  -  DELTA 
X(4)  =  MGT  -  DELTA 
X(ll)  =  CLT  -  3*DELTA 
X(39)  =  DELTA 
1030  RETURN 
END 

C - 

C - 


SUBROUTINE  TACHYHYDRITE 

C  Subroutine  for  CaMg2C16 . 12H20  equilibrium. 

REAL  IX, I, K 
CHARACTER* 50  TITLE 
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COMMON  /B/  X(60),GX(60),K(20:60),AH20,DEL(5), 

XIX(30) , TITLE, ITER,AX(60) , TEMP, I, PHI 

DOUBLE  PRECISION  KPRIME, CAT, MGT, CLT, X, DELTA 

DELTA=0 

CAT  =  X(3)  +  X(40) 

MGT  =  X(4)  +  2*X(40) 

CLT  =  X(ll)  +  6*X(40) 

KPRIME  =  K(40)/(GX(3)*GX(4)**2*GX(ll)**6*AH2O**12) 

IF  (CAT*MGT**2*CLT**6-KPRIME  .LT.  0)  GOTO  1050 
1040  X(3)  =  CAT 

X(ll)  =  CLT 
X(4)  =  MGT 
X(40)  =  0 
GOTO  1080 

1050  DECREM  =  MIN  (CAT,  MGT,  CLT-) 

DO  1070  KA  =  2,  5 

1060  DELTA  =  DELTA+DECREM* (DEL (KA) -1 ) 

IF  (DELTA  .LT.  0)  GOTO  1040 
X(3)  =  CAT  -  DELTA 

X(4)  =  MGT  -  2*DELTA 

X(ll)  =  CLT  -  6*DELTA 
IF  (KA  .EQ.  5)  CALL  PITZER 

KPRIME  =  K(40)/(GX(3)*GX(4)**2*GX(ll)**6*AH2O**12) 
IF  (X(3)*X(4)**2*X(11)**6-KPRIME  .GT.  0)  GOTO  1060 
DELTA  =  DELTA-DECREM* (DEL(KA)-l) 

1070  CONTINUE 

X(3)  =  CAT  -  DELTA 

X(4)  =  MGT  -  2*DELTA 

X(ll)  =  CLT  -  6*DELTA 
X(40)  =  DELTA 
1080  RETURN 
END 

C - 

C - 


SUBROUTINE  ARCANITE 

C  Subroutine  for  K2S04  equilibrium. 

REAL  IX, I, K 
CHARACTER*50  TITLE 

DOUBLE  PRECISION  KT, S04T, KPRIME, X, DELTA 

COMMON  /B/  X(60) ,GX(60) ,K(20:60) ,AH20,DEL(5) , 
xIX(30) , TITLE, ITER, AX(60) , TEMP, I, PHI 

DELTA=0 

KT  =  X(2)  +  2*X(45) 

S04T  =  X(12)  +  X(45) 

KPRIME  =  K(45) / (GX(2) **2*GX(12) ) 

IF  (KT**2*S04T-KPRIME  ,GT.  0)  GOTO  1100 
1090  X(2)  =  KT 

X(12)  =  S04T 
X(45)  =  0 
GOTO  1130 

1100  DECREM  =  MIN(KT,S04T) 

DO  1120  KA  =  2,  5 

1110  DELTA  =  DELTA+DECREM* (DEL (KA)-l) 

IF  (DELTA  .LT.  0)  GOTO  1090 
X(2)  =  KT  -  2*DELTA 
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X(12)  =  S04T  -  DELTA 
IF  (KA  .EQ.  5)  CALL  PITZER 
KPRIME  =  K(45)/(GX(2)**2*GX(12)) 

IF  (X(2)**2*X(12)-KPRIME  .GT.  0)  GOTO  1110 
DELTA  =  DELTA  -  DECREM* (DEL (KA) -1 ) 

1120  CONTINUE 

X(2)  =  KT  -  2*DELTA 
X(I2)  =  S04T  -  DELTA 
X(45)  =  DELTA 
1130  RETURN 
END 


C 

C 


SUBROUTINE  NA2S04 

C  Subroutine  for  Na2S04  equilibrium. 

REAL  MIRAB,IX,I,K 
CHARACTER*50  TITLE 

DOUBLE  PRECISION  NAT, S04T, KPRIME, X, DELTA 

COMMON  /B/  X(60) ,GX(60) ,K(20:60) ,AH20,DEL(5) , 
xIX(30) , TITLE, ITER, AX(60) , TEMP, I, PHI 

DELTA=0 

NAT  =  X(1)+2*X(41)+2*X(42) 

S04T  =  X(12)+X(41)+X(42) 

MIRAB  =  K(41) /AH2O**10 
TENAR  =  K(42) 

IF  (MIRAB  .GT.  TENAR)  THEN 

KPRIME  =  K(42) / (GX(1)**2*GX(12) ) 

ELSE 

KPRIME  =  K(41)/(GX(l)**2*GX(12)*AH2O**10) 

END  IF 

IF  (NAT**2*S04T-KPRIME  .GT.  0)  GOTO  1150 
1140  X(l)  =  NAT 

X(12)  =  S04T 
X(41)  =  0 
X(42)  =  0 
GOTO  1180 

1150  DECREM  =  MIN (NAT, S04T) 

DO  1170  KA  =  2,  5 

1160  DELTA  =  DELTA+DECREM* (DEL(KA) -1) 

IF  (DELTA  .LT.  0)  GOTO  1140 
X(l)  =  NAT  -  2*DELTA 
X(12)  =  S04T  -  DELTA 
IF  (KA  .EQ.  5)  CALL  PITZER 
MIRAB  =  K(41) /AH2O**10 
TENAR  =  K(42) 

IF  (MIRAB  .GT.  TENAR)  THEN 

KPRIME  =  K(42) / (GX(1)**2*GX(12) ) 

ELSE 

KPRIME  =  K(41)/(GX(1)**2*GX(12)*AH20**10) 

END  IF 

IF  (X (1) **2*X (12) -KPRIME  .GT.  0)  GOTO  1160 
DELTA  =  DELTA  -  DECREM* (DEL (KA) -1 ) 

1170  CONTINUE 

X(l)  =  NAT  -  2*DELTA 
X(12)  =  S04T  -  DELTA 
IF  (MIRAB  .LT.  TENAR)  THEN 
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X(41)  =  DELTA 
X(42)  =  0 

ELSE 

X(42)  =  DELTA 
X(41)  =  0 

END  IF 
1180  RETURN 
END 


C 

C 


SUBROUTINE  MGS04 

C  Subroutine  for  MgS04  equilibrium. 


REAL  IX, I, K 
CHARACTER* 50  TITLE 

DOUBLE  PRECISION  A, B, C, XA, MGT, S04T, KPRIME, KIONPR, X, DELTA 

COMMON  /B/  X(60) ,GX(60) ,K(20:60) ,AH20,DEL(5) , 

XlX(30) , TITLE, ITER, AX (60) , TEMP, I, PHI 

DELTA=0 

MGT  =  X(4)+X(43)+X(44)+X(23) 

S04T  =  X(12)+X(43)+X(44)+X(23) 

HEXAHY  =  K(43) /AH20**6 
EPSOM  =  K(44)  /  AH20**7 
•IF  (HEXAHY  .LT.  EPSOM)  THEN 

KPRIME  =  K(43)/(GX(4)*GX(12)*AH20**6) 

ELSE 

KPRIME  =  K(44)/(GX(4)*GX(12)*AH20**7) 

END  IF 

IF  (MGT*S04T-KPRIME  .GT.  O)  GOTO  1200 
1190  KIONPR  =  K(23)/(GX(4)*GX(12)) 

A  =  1 

B  =  -1* (MGT+S04T+KIONPR) 

C  =  MGT*S04T 

XA  =  (-B-SQRT(B**2-4*A*C))/(2*A) 

X(4)  =  MGT  -  XA 
X(12)  =  S04T  -  XA 
X(23)  =  XA 
X(43)  =  0 
X(44)  =  0 
GOTO  1230 

1200  DECREM  =  MIN(MGT,S04T) 

DO  1220  KA  =  2,  5 

1210  DELTA  =  DELTA+DECREM*(DEL(KA)-1) 

IF  (DELTA  .LT.  0)  GOTO  1190 
KIONPR  =  K(23) / (GX(4)*GX(12) ) 

X(23)  =  KPRIME/KIONPR 
X(4)  =  MGT  -  DELTA  -  X(23) 

X(12)  =  S04T  -  DELTA  -  X(23) 

IF  (KA  .EQ.  5)  CALL  PITZER 
HEXAHY  =  K(43) /AH20**6 
EPSOM  =  K(44) /AH20**7 
IF  (HEXAHY  .LT.  EPSOM)  THEN 

KPRIME  =  K(43)/ (GX(4)*GX(12)*AH20**6) 

ELSE 

KPRIME  =  K(44)/(GX(4)*GX(12)*AH20**7) 

END  IF 
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IF  (X(4)*X(12)-KPRIME  .GT.  0)  GOTO  1210 
DELTA  =  DELTA-DECREM* (DEL(KA) -1) 

1220  CONTINUE 

X(23)  =  KPRIME/KIONPR 
X(4)  =  MGT  -  DELTA  -  X(23) 

X(12)  =  S04T  -  DELTA  -  X(23) 

IF  (HEXAHY  .LT.  EPSOM)  THEN 
X(43)  =  DELTA 
X(44)  =  0 

ELSE 

X(44)  =  DELTA 
X(43)  =  0 

END  IF 
1230  RETURN 
END 


SUBROUTINE  PICROMERITE 

C  Subroutine  for  K2Mg (S04 ) 2 . 6H20  equilibrium. 

REAL  IX, I, K 
CHARACTER* 50  TITLE 

DOUBLE  PRECISION  A, B, C, XA, MGT, KT, S04T, KPRIME, KIONPR, X, DELTA 

COMMON  /B/  X(60),GX(60),K(20:60),AH20,DEL(5), 
xIX(30) , TITLE, ITER, AX (60) , TEMP, I, PHI 

PHASE  =  PHASE  +  1 
DELTA=0 

KPRIME  =  K(46)/(GX(4)*GX(2)**2*GX(12)**2*AH20**6) 

MGT  =  X(4)  +  X(46)  +  X(23) 

KT  =  X(2)  +  2*X(46) 

S04T  =  X(12)  +2*X(46)  +  X(23) 

IF  (MGT*KT**2*S04T**2-KPRIME  .GT.  0)  GOTO  1250 
1240  KIONPR  =  K(23)/(GX(4)*GX(12)) 

A  =  1 

B  =  -1* (MGT+S04T+KIONPR) 

C  =  MGT*S04T 

XA  =  (-B-SQRT(B**2-4*A*C) )/ (2*A) 

X(4)  =  MGT  -  XA 
X(12)  =  S04T  -  XA 
X(23)  =  XA 
X(2)  =  KT 
X(46)  =  0 
GOTO  1280 

1250  DECREM  =  MIN (MGT, S04T, KT) 

DO  1270  KA  =  2,  5 

1260  DELTA  =  DELTA+DECREM* (DEL (KA) -1) 

IF  (DELTA  .LT.  0)  GOTO  1240 
KIONPR  =  K(23) / (GX(4) *GX(12) ) 

X(23)  =  KPRIME/KIONPR 
X(4)  =  MGT  -  DELTA  -  X(23) 

X(12)  =  S04T  -  2*DELTA  -  X(23) 

X(2)  =  KT  -  2*DELTA 
IF  (KA  .EQ.  5)  CALL  PITZER 

KPRIME  =  K(46)/(GX(4)*GX{2)**2*GX(12)**2*AH20**6) 

IF  (X(4)*X(2)**2*X(12)**2-KPRIME  .GT.  0)  GOTO  1260 
DELTA  =  DELTA-DECREM* (DEL (KA)-l) 

1270  CONTINUE 
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X(23)  =  KPRIME/KIONPR 
X(4)  =  MGT  -  DELTA  -  X(23) 

X(12)  =  S04T  -  2*DELTA  -  X(23) 

X(2)  =  KT  -  2*DELTA 
X(46)  =  DELTA 
1280  RETURN 
END 

C— - 

C - 

SUBROUTINE  CAS04 

C  Subroutine  for  CaS04  equilibrium. 

REAL  IX, I, K 
CHARACTER*50  TITLE 

DOUBLE  PRECISION  A, B, C, XA, CAT, S04T, KIONPR, IXA, X 

COMMON  /B/  X(60) ,GX(60) ,K(20:60) ,AH20,DEL(5) , 
xIX(30) , TITLE, ITER, AX (60) , TEMP, I,  PHI 

1290  CALL  PITZER 

KIONPR  =  K(22) / (GX(3) *GX(12) ) 

CAT  =  X(3)  +  X(22) 

S04T  =  X(12)  +  X(22) 

IXA  =  X(22) 

A  =  1 

B  =  -1*(CAT  +  S04T  +  KIONPR) 

C  =  CAT*S04T 

XA  =  (-B-SQRT(B**2-4*A*C) ) / (2*A) 

X(3)  =  CAT  -  XA 
X(12)  =  S04T  -  XA 
X(22)  =  XA 

IF  ( (ABS (IXA-XA) /XA) *100  .GT.  1)  GOTO  1290 

RETURN 

END 

C - 
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